dolfin team mailing list archive

[Bartosz Sawicki <sawickib@xxxxxxxxxxxxx>]

Anders Logg:
> On Mon, Mar 31, 2008 at 02:42:14PM +0200, Kristen Kaasbjerg wrote:
>> Anders Logg wrote:
>>> On Mon, Mar 31, 2008 at 02:26:45PM +0200, Kristen Kaasbjerg wrote:
>>>   
>>>> Hi,
>>>>
>>>> I would like to plot a function that is the difference between a 
>>>> DiscreteFunction and a UserFunction. Can this be done in any smart ways ?
>>>>
>>>> The only way I can think of is to create a new user defined function 
>>>> that contains the two Functions above, which does the subtraction in the 
>>>> eval function, and then save this function to vtk format.
>>>> This, however, will call the eval method of the DiscreteFunction which 
>>>> includes cell searching etc. Just sounds a bit cumbersome to me when the 
>>>> mesh contains say 500.000 elements.
>>>>
>>>> Kristen
>>>>     
>>> I would recommend to project the difference to an appropriate finite
>>> element space and then plot the projection.
>>>
>>>   
>> Ok, are there any examples around that illustrate how to do this ?
>>
>> Kristen
> 
> v = TestFunction(...)
> u = TrialFunction(...)
> f = Function(...)
> g = Function(...)
> 
> a = v*u*dx
> L = v*(f - g)*dx
> 

I had similar problem, where vector elements was used:

ve = VectorElement("Discontinuous Lagrange", "tetrahedron", 0)
v = TestFunction(ve)
err = TrialFunction(ve)
anal = Function(ve)
sol = Function(ve)
a = dot(v,err)*dx
L = dot(v, (anal-sol) )*dx

This example works, but form language is sometimes not enough.
I've failed for following form:

e = FiniteElement("Discontinuous Lagrange", "tetrahedron", 0)
ve = VectorElement("Discontinuous Lagrange", "tetrahedron", 0)
v = TestFunction(e)
err = TrialFunction(e)
anal = Function(e)
sol = Function(ve)
a = v*err*dx
L = v*modulus(anal - sqrt(sol[0]*sol[0]+sol[1]*sol[1]+sol[2]*sol[2]))*dx

*** Error at w1_a8w1_a9 | va8[0]*va9[0] + w1_a12w1_a13 | va12[1]*va13[1]
+ w1_a16w1_a17 | va16[2]*va17[2]
*** Cannot take square root of sum.

Finally I decided to write simple loop over DiscreteFunction vector
elements. Works perfectly even for 1M elements meshes :)

cheers
BArtosz