ffc team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Tuesday 25 November 2008 20:51:58 Anders Logg wrote:
> On Tue, Nov 25, 2008 at 03:56:19PM +0100, Alessio Quaglino wrote:
> > There is a function in FFC which seems very useful for my aim but I don't
> > understand how to use it.  One part of my form is the following:
> >
> > n = FacetNormal("triangle")
> > dsD = Integral("exterior facet", 0)
> > dsN = Integral("exterior facet", 1)
> > dsF = Integral("exterior facet", 2)
> >
> > a = ( dot(psi,b) - div(psi)*w + phi*div(b) )*dx + ( dot(b,n)*psi[0] +
> > dot(b,n)*psi[1] )*dsF L = phi*f*dx - dot(psi,g)*dsN
> >
> > What I want to do is to assemble the matrix splitting the boundary in 3
> > parts: dsD (no terms), dsN (1 term in a), dsF (1 term in L). However I
> > don't understand how to pass to the form the subdomains I defined in
> > dolfin. For example I have:
> >
> > class DirichletBoundary : public SubDomain
> > {
> >   bool inside(const double* x, bool on_boundary) const
> >   {
> >     return x[0] < DOLFIN_EPS;
> >   }
> > };
> >
> > which I would like to associate to dsD. Is there a way to do that?
> >
> > Thank you,
> > Alessio
>
> Take a look at demo/pde/lift-drag/cpp/main.cpp in DOLFIN.

I do not know if that is a good example for you case.

In that example a SubDomain marking the subdomain that shall be assembled is 
passed to the assemble function. This wont work in your case. As the passed 
SubDomain will be used to create a MeshFunciton over both facet and cells. I 
assume that your SubDomain is created to define inside if you are on a 
special boundary, and the resulting MeshFunction defining the cells, will 
then be defined as "outside". Please correct me here if I am totaly of 
Anders.

For your need I would create the actuall meshfunctions, defining the 
subdomains. I assume you have defined three SubDomains called sd0, sd1, sd2, 
all overloading the inside function that defines the your subdomain, as you 
have above. Also see demo/mesh/subdomains.

    ....
    // Creating a dummy cell_domain MeshFunction, maybee an zero pointer
    // would do the job?
    MeshFunction<uint> cell_domains(mesh, mesh.topology().dim());
    cell_domains = 0;
    
    MeshFunction<uint> facet_domains(mesh, mesh.topology().dim() - 1);
    // First mark all subdomains with an "outside" marker, defining
    // the rest of the domain.
    facet_domains = 3;
    sd0.mark(facet_domains,0);
    sd1.mark(facet_domains,1);
    sd2.mark(facet_domains,2);
    
    Matrix A;
    assemble(A,form,mesh,cell_domains,facet_domains,facet_domains);   

I am not a native c++ dolfin user so I might have got some of the syntax wrong 
above, but I hope you get the point. 

Note that you send the facet_domains as sub_domain markes for the interiour 
facet integrals too, see signature of assemble() in assemble.h, but that wont 
be used as no such integral is defined in your form.



Johan