dolfin team mailing list archive

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

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
 
> -- 
> Anders
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>