dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Thursday 12 June 2008 22:05:43 Anders Logg wrote:
> On Wed, Jun 11, 2008 at 08:27:14AM +0200, Johan Hake wrote:
> > On Tuesday 10 June 2008 14:52:52 Anders Logg wrote:
> > > On Tue, Jun 10, 2008 at 02:28:51PM +0200, Johan Hake wrote:
> > > > On Tuesday 10 June 2008 13:23:45 Anders Logg wrote:
> > > > > On Tue, Jun 10, 2008 at 09:55:51AM +0200, Martin Sandve Alnæs wrote:
> > > > > > 2008/6/9 DOLFIN <dolfin@xxxxxxxxxx>:
> > > > > > > One or more new changesets pushed to the primary dolfin
> > > > > > > repository. A short summary of the last three changesets is
> > > > > > > included below.
> > > > > > >
> > > > > > > changeset:   4282:5f6759a7460d81435abf4d67911667e3e628e164
> > > > > > > tag:         tip
> > > > > > > user:        Anders Logg <logg@xxxxxxxxx>
> > > > > > > date:        Mon Jun 09 23:06:51 2008 +0200
> > > > > > > files:       dolfin/fem/Assembler.cpp dolfin/fem/UFC.cpp
> > > > > > > dolfin/fem/UFC.h description:
> > > > > > > Reset local tensor A^K before call to tabulate_tensor()
> > > > > >
> > > > > > Why? This is the job of tabulate_tensor, which can do it much
> > > > > > more efficiently?
> > > > >
> > > > > I added this for convenience. It is used when assembling forms over
> > > > > subdomains and the form compiler needs to generate code for zero
> > > > > contributions from other domains, for example when integrating
> > > > >
> > > > >   M = f*ds(5)
> > > > >
> > > > > Then the compiler needs to add zero integrals for ds([0-4]).
> > > > > Then this code just needs to be empty.
> > > >
> > > > It would be great to be able to integrate the above mentioned form.
> > > > It took me some time to realise that you could only assemble
> > > > subdomains that was contiguously numbered from 0 and up.
> > >
> > > That's not needed anymore. It should be possible now to integrate the
> > > form above without adding any additional terms.
> >
> > Nice!
> >
> > > > But is this implementation the most logical? What if I think it is
> > > > natural to mark a subdomain with "100"? Then 100 empty integrals
> > > > would be created? Is it possible to store the integrals in a map,
> > > > with markers as keys. Or store the integrals in a contiguous array,
> > > > and then map these to markers?
> > >
> > > Yes, this can be done. It's related to building a mapping from
> > > subdomains to cells that's been discussed before. Now we only have the
> > > inverse, the mapping from cells to subdomains. I don't know when this
> > > will happen. It's something I'd like to have but it's not critical.
> >
> > OK
> >
> > > > This would add flexibility for assigning a form integral an arbitrary
> > > > number, but it also introduce a map lookup overhead that would slow
> > > > down the assembly. The latter could be solved with an inverse
> > > > MeshFunction, which map values to cells/facets.
> > >
> > > Don't think it needs to slow down assembly.
> >
> > OK
> >
> > > > It would also be nice to separate a facet MeshFunction in one only
> > > > for interior facets and one only for the exterior ones. This would
> > > > fit nicely in to the above mentioned procedure. I do not know how to
> > > > do this nicely though...
> > >
> > > Yes.
> > >
> > > > A note related to this issue:
> > > > A typical small mesh I generate have 23K vertices, and 86K
> > > > tetrahedrons, it takes ~0.3 s to load. The boundary mesh have 35K
> > > > facets but the whole mesh have 190K facets, including the interior
> > > > facets. This takes ~5s to load, almost 20 times more time than the
> > > > mesh, and it loads information about 160K facets that is not needed.
> > >
> > > I guess you need this for boundary conditions? If that is the case,
> > > then the internal facets never need to be generated if you use the new
> > > mesh file format (which VTK uses) where all the boundary information
> > > is included in the mesh file.
> >
> > ?
> >
> > Maybee I was unclear? I need a MeshFunction to define my boundary
> > subdomains, to send to assemble. Do you say that I can store this
> > subdomain meshfunction in the new mesh format, by only storing the
> > information of the exterior facets? I thought a Meshfunction was tied to
> > the same entities in the mesh, and needed to be of the same length, i.e.,
> > a facet meshfunction need to include information of all facets in the
> > mesh.
>
> Yes, that's true, but in the "new" (extended) mesh format, we store
> Arrays, not MeshFunctions. These don't have to have the same length.
>
> Take a look at data/meshes/aneurysm.xml.gz.
>
> Only DirichletBC can currently handle this data. Assembler still
> relies on a MeshFunction for the subdomains.

Ok, this MeshData can take an arbitrary number of MeshFunctions and Arrays. 
But three arrays with name "boundary facet cells", "boundary facet numbers" 
and "boundary indicators" are treated specially in the function DirichleBC().

I do not really understand what these name stand for. Shouldn't it be enough 
to use two arrays? One for the global facet number of the facets laying on 
the boundary, and one for the markers?

I still can't figure out why a mesh takes ~20 times logner to load? First I 
thought that a MeshFunction stores more information, but 20 times, and the 
file is actually much smaller than the mesh file. Could the implementation of 
XMLMeshFunction be to blame? 

A MeshFunction for a slightly larger mesh, than the above mentioned, takes 
~1min while the mesh takes 2.5 s to load... 

Try for your self:
-------------------------------------------------
wget http://folk.uio.no/hake/mesh/\
{face_markers.xml.gz,mesh_test.py,mesh.xml.gz}
python mesh_test.py
-------------------------------------------------

Johan