dolfin team mailing list archive
-
dolfin team
-
Mailing list archive
-
Message #20234
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