dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Dec 11, 2007 at 09:58:37AM +0100, Kristian Oelgaard wrote:
> Quoting Anders Logg <logg@xxxxxxxxx>:
> 
> > On Tue, Dec 04, 2007 at 07:55:31PM +0100, Kristian Oelgaard wrote:
> > 
> > > Hi,
> > > 
> > > I'm having a few difficulties with 1D elements in DOLFIN, more specifically
> > with
> > > the assembly over exterior integrals and when applying boundary
> > conditions.
> > > 
> > > In Assembler::assembleExteriorFacets(), a boundary mesh is being created.
> > This
> > > fails when the mesh consists of intervals. It does so in CellType where
> > only
> > > Interval, Triangle and Tetraheron are available.
> > > 
> > > I successfully implemented a special case in
> > Assembler::assembleExteriorFacets()
> > > if (mesh.topology().dim() == 1){//handle 1D}, but is this what we want and
> > could
> > > it be done more naturally?
> > 
> > I think we can add a new cell type for 0-D cells. The natural thing
> > would be to call it Point, but it's already used in DOLFIN for
> > something else. Maybe we should rename Triangle to TriangleCell etc
> > since a Triangle is not really a Triangle. It's just a collection of
> > functions on Triangle cells. Then we can have PointCell, IntervalCell,
> > TriangleCell, TetrahedronCell.
> 
> Sounds like an OK solution to me. Should I take a look at this?

That would be excellent.

> > > The second problem (BCs) is due to how VertexIterator is defined.
> > > Using the 'topological' (default) or 'geometrical' method the following
> > simple
> > > SubDomain will fail:
> > > 
> > > class DirichletBoundary : public SubDomain
> > > {
> > >   bool inside(const real* x, bool on_boundary) const
> > >   {
> > >     return std::abs(x[0]) < DOLFIN_EPS && on_boundary;
> > >   }
> > > };
> > > 
> > > When a facet with index 0 (this is the facet I want) is used to generate a
> > > VertexIterator it will look at the connectivity and return the vertex with
> > index
> > > 1 on a standard mesh. A VertexIterator generated from the facet with index
> > 1
> > > will iterate over vertices 0 and 2. This destroys the logic of the
> > geometrical
> > > check in the inside() function.
> > 
> > I don't understand this. Which mesh is this on?
> 
> On an interval mesh.
>   0     1     2
>   *-----*-----*
> 
> SubDomain.mark()
> 
> will loop all facets, in this case 0,1,2, and check if a facet is a boundary
> facet i.e. if it has only 1 neighboring cell. This test is OK. It will then
> create a VertexIterator to check if all vertices are on the boundary. This
> iterator will be created from the 0 - 0 connectivity so in this case:
> 
> 0: 1
> 1: 0 2
> 2: 1
> 
> facet 0, on_boundary = true: check if vertex 1 is on boundary, false   NOT OK
> facet 2, on_boundary = true: check if vertex 1 is on boundary, false   NOT OK
> 
> So in this case no boundaries will be marked.

ok, I see. Just modify it like this:

  // Dimension of facet > 0, check incident vertices
  if (entity->dim() > 0)
  {
    for (VertexIterator vertex(*entity); !vertex.end(); ++vertex)
    {
      simple_array<real> x(mesh.geometry().dim(), vertex->x());
      if (!inside(x, on_boundary))
      {
        all_vertices_inside = false;
        break;
      }
    }
  }
  // Dimension of facet = 0, so just check the vertex itself
  else
  {
    simple_array<real> x(mesh.geometry().dim(), mesh.geometry().x(entity->index()));
    if (!inside(x, on_boundary))
    {
      all_vertices_inside = false;
      break;
    }
  }

-- 
Anders