dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

There are two things that need to happen for boundary conditions
(or boundary integrals).

1. FFC needs to generate the code for computation of element tensors
on boundaries. When FFC compiles a form, all terms containing '*dx'
are compiled to generate the function

  void BilinearForm::eval(real block[], const AffineMap& map) const;

FFC needs to be updated to compile all terms containing '*ds' into the
function

  void BilinearForm::eval(real block[], const AffineMap& map, uint segment) const;

where the extra argument segment tells the eval() function on which
part of the boundary of the current cell (triangle or tetrahedron) the
element tensor (element stiffness matrix) is to be computed.

2. DOLFIN needs to loop over the boundary and call this function on
each cell located at the boundary during assembly. This should happen
in the assembly functions in the class FEM (src/kernel/fem/FEM.cpp) in
DOLFIN. I suggest adding a new function FEM::assembleBoundary() to
this class that can be called from the assembly functions. We also
need two functions assembleBoundary_2D() and assembleBoundary_3D(),
since we need to do things differently for triangles and tetrahedra.

The assembleBoundary() functions will look similar to the applyBC()
functions. For triangles, we iterate over the edges of the boundary,
then for each edge pick the triangle that the edge is part of (there
is only one such triangle since the edge is on the boundary). We also
need to check which edge of the triangle we are integrating over. The
edge can be edge number 0, 1 or 2 of the triangle (this is the
argument segment to eval()). There is already a class AffineMap in
DOLFIN that maps from the reference triangle, but I think we need to
either extend this map or implement a new one that takes into account
which edge we are integrating over.

The code would look something like

for (EdgeIterator edge(boundary); !edge.end(); ++edge)
{
    // Get cell containing the edge (pick first, should only be one)                                     dolfin_assert(edge->noCellNeighbors() == 1);
    const Cell& cell = edge->cell(0);

    // Compute the local edge number (0, 1 or 2)
    segment = ...

    // Update affine map (maybe change to another map for boundaries)
    map.update(cell);

    // Compute map from local to global degrees of freedom
    test_element.dofmap(test_dofs, *cell, mesh);
    trial_element.dofmap(trial_dofs, *cell, mesh);

    // Compute element matrix
    a.eval(block, map, segment);

    // Add element matrix to global matrix
    A.add(block, test_dofs, m, trial_dofs, n);
}

To get going, you could take a simple problem and implement the
function BilinearForm::eval() for boundaries by hand and then
implement the boundary assembly in FEM.cpp. Then I can implement
the generation of the eval() function in FFC as soon as I find the
time.

Tell me what you think and get back with more questions.

/Anders

On Thu, Jun 30, 2005 at 02:48:06PM +0200, Karin Kraft wrote:
> Hello,
> 
> We started to look at the boundary conditions yesterday. Do you have any
> plans on how to do it? I have started to implement a mapping that takes
> a triangle in 3D to a unit triangle in 2D but I am not sure of the
> consequences for the rest of the code so I am happy for some
> suggestions.
> 
> /Karin
> 
> On Wed, 2005-06-08 at 08:53 -0500, Anders Logg wrote:
> > Hi Karin!
> > 
> > The only thing we have right now is Dirichlet boundary conditions or
> > homogeneous Neumann boundary conditions. To do anything else, we need
> > to be able to compute integrals over the boundary and there is nobut I am not sure of the consequences for the rest of th
> > support for this in yet in FFC or DOLFIN.
> > 
> > It won't be very difficult to add, once I find time to do it. If
> > someone else wants to try, you are more than welcome. It would involve
> > modifiying both FFC and DOLFIN. Tell me if you want some detailed
> > pointers.
> > 
> > In the meantime, it might be possible to cheat by approximating a
> > boundary integral with an integral over the domain, using a function
> > that is zero everywhere except on elements close to the boundary,
> > where the function takes the value ds/dx (length of boundary edge
> > divided by area of triangle for a triangular mesh).
> > 
> > /Anders
> > 
> > On Wed, Jun 08, 2005 at 03:06:29PM +0200, Karin Kraft wrote:
> > 
> > > Hello,
> > > 
> > > I am wondering how the support for boundary terms is in Dolfin. I would
> > > like to use absorbing boundary conditions. Johan Jansson said that at
> > > the moment it is not working. Is somebody working on it? Otherwise, how
> > > difficult is it to add? 
> > > 
> > > Regards,
> > > Karin Kraft
> > > 
> > > 
> > > _______________________________________________
> > > DOLFIN-dev mailing list
> > > DOLFIN-dev@xxxxxxxxxx
> > > http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> > > 
> > 
> 
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> 

-- 
Anders Logg
Research Assistant Professor
Toyota Technological Institute at Chicago
http://www.tti-c.org/logg/