dolfin team mailing list archive

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

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


 
> -- 
> Anders
> 
> 
> > However, using the 'pointwise' method the above SubDomain works as
> expected.
> > 
> > Any thoughts on where these problems can be tackled most efficiently?
> > 
> > Kristian
> > 
> > 
> > 
> > 
> > 
> > 
> > _______________________________________________
> > DOLFIN-dev mailing list
> > DOLFIN-dev@xxxxxxxxxx
> > http://www.fenics.org/mailman/listinfo/dolfin-dev
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev
>