dolfin team mailing list archive

[Andre Massing <massing@xxxxxxxxx>]

Den 22.11.10 18.39, skrev Anders Logg:
> On Mon, Nov 22, 2010 at 05:45:07PM +0100, Andre Massing wrote:
> > Den 22.11.10 17.40, skrev Johan Hake:
> > > On Monday November 22 2010 07:36:25 Andre Massing wrote:
> > > > Hi!
> > > >
> > > > I am not entirely sure, whether there was a similar discussion about that,
> > > >    but I just recognized that using a ufl file like
> > > >
> > > > element = FiniteElement("Lagrange","tetrahedron", 1)
> > > > w = Coefficient(element)
> > > > M = w*w*dx(0)
> > > >
> > > > and a meshfunction with 0, 1, 2 values
> > > > results into the assembly of a function on *entire* domain and not only
> > > > cells that a marked with 0.
> > > > Can anybody confirm this behaviour?
> Take a look in Assembler.cpp in the function assemble_cells. It should
> be quite easy to find what goes wrong with some suitable print
> statements. Look for the if-statements starting with
>
>    if (domain&&  domains->size()>  0)
I cannot really find the culprit code snippets. Either I am completely 
blind right now, or it is rather strange...
I have an example where I use a subdomain directly (meaning a SubDomain 
class instance) in the assemble
routine and then indirectly, via marking a meshfunction via subdomain 
(which is the what the assemble function
does it internally). This gives different result and I cannot see why.
Code is attached...

--
Andre


#Form to compile L^2.

element = FiniteElement("Lagrange","tetrahedron", 1)
w = Coefficient(element)
M = w*w*dx(0)

element = FiniteElement("Lagrange", "tetrahedron", 1)

v = TestFunction(element)
u = TrialFunction(element)
f = Coefficient(element)

a = inner(grad(v), grad(u))*dx
L = v*f*dx

#include "L2_norm.h"
#include "Poisson.h"
#include <dolfin.h>

using namespace dolfin;

class HalfSpace : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
    return x[0] < 0.5 + DOLFIN_EPS;
  }
};


int main()
{
    UnitCube mesh0(5,5,5);

    HalfSpace half_space;
    MeshFunction<uint> half_space_marker(mesh0, mesh0.topology().dim(),1);
    half_space.mark(half_space_marker,0);
    info(half_space_marker.str(true));

    Poisson::FunctionSpace V0(mesh0);
    Function one(V0);
    one.interpolate(Constant(1));

    L2_norm::Functional one_norm(mesh0,one);

    double norm;
    norm = assemble(one_norm,half_space);
    info("Integral of 1 is %e:", norm);

    norm = assemble(one_norm,&half_space_marker);
    info("Integral of 1 is %e:", norm);


    return 0;
}