← Back to team overview

dolfin team mailing list archive

Re: assembly of whole domain if domain marker = 0 is used

 

The problem is you are calling the wrong assemble function. Look in
assemble.h. You are calling this function:

  double assemble(const Form& a,
                  bool reset_sparsity=true,
                  bool add_values=false);

and the pointer of the marker is interpreted as a boolean.

You want this function:

  double assemble(const Form& a,
                  const MeshFunction<uint>* cell_domains,
                  const MeshFunction<uint>* exterior_facet_domains,
                  const MeshFunction<uint>* interior_facet_domains,
                  bool reset_sparsity=true,
                  bool add_values=false);

Change to

  norm = assemble(one_norm, &half_space_marker, 0, 0);

and it should work.

--
Anders



On Mon, Nov 22, 2010 at 08:36:19PM +0100, Andre Massing wrote:
> 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...
>

> #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;
> }

> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp




Follow ups

References