dolfin team mailing list archive

[Kristian Oelgaard <k.b.oelgaard@xxxxxxxxxx>]

Quoting 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.

On second thought, the current definition of VertexIterator makes most sense. If
one already has the facet/vertex of interest there's no need to create an
iterator that essentially would be the same entity.

I'll fix it.

Kristian
 
> -- 
> Anders
>