dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Dec 11, 2007 at 10:32:56AM +0100, Kristian Oelgaard wrote:
> Quoting 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.
> 
> OK
>  
> > > > > 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;
> >     }
> >   }
> 
> Sure, that would fix this problem, but what about all other cases where a
> VertexIterator is created from a MeshEntity and this entity happens to be a
> Vertex itself? Would it make more sense to modify the VertexIterator?
> 
> Kristian

The current definition is that the incident vertices of a vertex are
those vertices that share a common cell (and also edge and face for
simplices), not including the vertex itself. I think that makes most
sense but I'm open to suggestions/comments.

-- 
Anders