← Back to team overview

dolfin team mailing list archive

Re: [HG DOLFIN] Reset local tensor A^K?before?call?to?tabulate_tensor()

 

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().)

-- 
Anders

Attachment: signature.asc
Description: Digital signature


Follow ups

References