dolfin team mailing list archive

[Karin Kraft <karin.kraft@xxxxxxxxxxx>]

Hello!

We have made some progress. We have made a simple implementation of the
general boundary conditions. As a start we have assumed homogeneous
dirichlet with a constant penalty factor and this works reasonably
well. 

Some issues we have discovered:
1. We need to compute the local edge numbers. At the moment we do this
by comparing to all the edges in the cell like so:

  // Iterate over all edges in the boundary
  for (EdgeIterator edge(boundary); !edge.end(); ++edge)
  {
      ...
      for(int i = 0; i < cell.noEdges(); i++)
      {
	  if(cell.edgeID(i) == edge->id())
	  {
	      segment = i;
	  }
      }


However, this is inefficient. Perhaps the mesh should be extended to
store this information?

2. Our eval looks something like this and we want to avoid the
conditional evaluation. Any ideas?
 
 void eval(real block[], const AffineMap& map, uint boundary) const
  {
    // Compute geometry tensors
    real G0_ = map.scale*1e6;

    // Compute element tensor

    if( boundary == 0 ) {
	block[0] = 0.0;
	block[1] = 0.0;
	block[2] = 0.0;
	
	block[3] = 0.0;
	block[4] = 3.333333333333318e-01*G0_;
	block[5] = 3.333333333333318e-01*G0_;

	block[6] = 0.0;
	block[7] = 3.333333333333318e-01*G0_;
	block[8] = 3.333333333333318e-01*G0_;
    }

Regards,
Karin and Johan

On Thu, 2005-06-30 at 11:57 -0500, Anders Logg wrote:
> 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
> > 
>