dolfin team mailing list archive

[Johan Hake <hake@xxxxxxxxx>]

On Sunday 15 June 2008 00:22:33 Anders Logg wrote:
> On Fri, Jun 13, 2008 at 12:00:53PM +0200, Johan Hake wrote:
> > On Friday 13 June 2008 10:57:00 Anders Logg wrote:
> > > On Fri, Jun 13, 2008 at 10:43:12AM +0200, Johan Hake wrote:
> > > > On Friday 13 June 2008 10:17:20 Anders Logg wrote:
> > > > > On Fri, Jun 13, 2008 at 10:08:31AM +0200, Johan Hake wrote:
> > > > > > 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().
> > > > >
> > > > > Correct.
> > > > >
> > > > > > 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?
> > > > >
> > > > > No, it's not enough. It would be enough if everyone numbered the
> > > > > facets in the same way as DOLFIN, but that is not the case, in
> > > > > particular since the way DOLFIN numbers the facets depends on the
> > > > > particular algorithm implemented in DOLFIN for generating the
> > > > > facets for a mesh. So we need to store the markers based on the
> > > > > *cells* rather than the *facets*. So we basically store a list of
> > > > > the cells adjacent to the facets on the boundary, the local
> > > > > numbering of the facets relative to those cells, and the markers
> > > > > for the facets.
> > > >
> > > > OK!
> > > >
> > > > > > 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?
> > > > >
> > > > > I think the problem is that when you load markers on the facets of
> > > > > a mesh, then DOLFIN first needs to compute the facets and that may
> > > > > take some time. (Only the vertices and cells are present in the XML
> > > > > file.)
> > > >
> > > > Jepp, you are correct here!
> > > >
> > > > mesh.init(2) took 55 s and loading the mesh function after this took
> > > > 2s... :P
> > >
> > > Good. (But actually not that good. Looks like we might need to do
> > > something to improve the speed of mesh.init().)
> >
> > Yes it takes quite a while...
> >
> > I timed the TopologyComputation::computeEntities() function:
> > Start initializing mesh:
> >   ComputeConnectivity: 25.60s
> >   Initializing local arrays: 0.00s
> >   Counting entities: 5.94s
> >   Initialising entities and connections: 0.04s
> >   Adding entities: 23.60s
> >   Deleting data: 0.00s
> > Initializing mesh took 55.35s
> >
> > One could probably get a better view with a propper profiling. But I just
> > made a quick timing of the different parts of the function. Looks like
> > most of the time is spent in initalising the connectivities and creating
> > and add the actuall entities.
>
> I haven't made any attempt to optimize any part of the implementation
> in TopologyComputation. I don't think anything can be done about the
> complexity but one could probably do some clever tricks to speed it
> up.
>
> Anyway, I don't have time to dig into this now. Feel free to take a
> look or else post it as a feature request in Bugzilla.

I filed a feature request at Bugzilla.

Johan