Subject: [FFC-dev] subdomains

From: Alessio Quaglino <[email-masked]>
Date: Tue, 25 Nov 2008 15:56:19 +0100

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

---

From: Anders Logg <[email-masked]>
Date: Tue, 25 Nov 2008 20:51:58 +0100

> which I would like to associate to dsD. Is there a way to do that?

Take a look at demo/pde/lift-drag/cpp/main.cpp in DOLFIN.

-- 
Anders

---

From: Johan Hake <[email-masked]>
Date: Wed, 26 Nov 2008 10:07:32 +0100

On Tuesday 25 November 2008 20:51:58 Anders Logg 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:
> >
> > 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:
> >

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

---

From: Anders Logg <[email-masked]>
Date: Wed, 26 Nov 2008 12:03:17 +0100

> then be defined as "outside". Please correct me here if I am totaly of 

Correct. You will need to define a MeshFunction.

> be used as no such integral is defined in your form.

Take a look at demo/mesh/subdomains/cpp/main.cpp which demonstrates
how to use SubDomain to mark sub domains in a MeshFunction.

Note that the order in which you mark sub domains is important (if
they are not disjoint).

-- 
Anders