> Date: Wed, 5 Jun 2013 08:31:50 +0200
> From: [email protected]
> To: [email protected]
> CC: [email protected]
> Subject: Re: [FEniCS] FEM LIBRARY WITH OCTAVE
>
> On 06/04/2013 11:25 PM, Marco Vassallo wrote:
> >> Date: Thu, 30 May 2013 09:47:27 +0200
> >> From: [email protected]
> >> To: [email protected]
> >> CC: [email protected]
> >> Subject: Re: [FEniCS] FEM LIBRARY WITH OCTAVE
> >>
> >> On Thu, May 30, 2013 at 09:42:12AM +0200, Martin Sandve Alnæs wrote:
> >> > On 30 May 2013 09:34, Johan Hake <[1][email protected]> wrote:
> >> >
> >> > Hi Marco!
> >> > DOLFIN use SWIG to create Python wrappers. SWIG provides Octave
> >> > wrappers
> >> > too. Are you going to use those?
> >> > [2]http://www.swig.org/Doc2.0/SWIGDocumentation.html#Octave
> >> > SWIG use a flat namespace for the generated wrapper code, so you
> >> > would
> >> > still have nameclashes for the Array class. However it will
> >> > internally
> >> > keep all original C++ type information. SWIG also allows you to
> >> > rename
> >> > any C++ entity, classes, functions, variables, methods, aso, to
> >> > whatever
> >> > in the target language. You can then just tell SWIG to rename the
> >> > dolfin::Array to, let say, dolfin::ArrayDolfin.
> >> > Another challenge you might consider before you invest too much is
> >> > to
> >> > figure out how you should integrate the code generation part of
> >> > FEniCS.
> >> > For now that is handled by the Python packages of FFC, FIAT, UFL,
> >> > and
> >> > parts of UFC. Because these are all in Python we have been able to
> >> > integrate these seamlessly (but not without bumps) with the Python
> >> > interface of DOLFIN.
> >> >
> >> >
> >> > Doing this within octave will be challenging.
> >> >
> >> > To say the least... I'd say it's close to impossible to get the same
> >> > level of integration with just in time compilation as we have in
> >> > python.
> >> > However if it's possible to wrap pregenerated forms (i.e. header files
> >> > generated by running ffc on the commandline) it could be useful to
> >> > assemble matrices directly in octave for analysing eigenvalues etc. In
> >> > some older internal code here at Simula we had the concept of a "matrix
> >> > factory" which could assemble common matrices from a number of
> >> > pregenerated form files, e.g. mass, stiffness, convection etc.
> >>
> >> A matrix factory approach sounds like a very good fit for Octave. I
> >> suggest using that approach and implementing an extensive list of
> >> "finite element matrices" for various degrees of polynomials and
> >> possibly different element families.
> >>
> > Following your advice, as a first step for the matrix factory I'm creating
> > a function for the conversion of .xml mesh to/from (p,e,t) format.
> > 1) In the conversion from (p,e,t) to xml I'm doing something very
> > simple like
> >
> > for(uint i = 0; i < num_vertex; ++i)
> > editor.add_vertex(i, p(0,i), p(1,i));
> >
> > for(uint i = 0; i < num_cells; ++i)
> > editor.add_cell(i, t(0,i)-1, t(1,i)-1, t(2,i)-1);
> >
> > which works fine. But I have some problem when i try to use the information
> > contained in the "e" matrix about the number of the geometrical border
> > containing the side edge.
>
> What is contained in e? Is it edge information? A DOLFIN mesh does only
> use the point-triangle (vertex-cell) connectivity. Edge connectivities
> can be generated from these connectivities.
>
> > I would like to use something like
> >
> > FacetFunction<uint> boundary_markers(mesh);
> >
> > as it is showed in the book.
> > I would like to assign the information about the geometrical border
> > while I'm editing the mesh instead of iterating over all the facets of
> > the mesh and then recover the information with another loop on the
> > "e" matrix looking for the corresponding Facets. (And even in this case,
> > how can I get extract information of the indexes related to the facets?)
>
> You can't. A MeshFunction (which FacetFunction is) can only be
> constructed once you have made your Mesh and generated the
> connectivities mentioned above.
>
> > Is something like this possible??
> > for(uint i = 0; i < num_side_edges; ++i){
> > editor.add_facets(i, from vertex a(i) to b(i));
> > boundary_markers[i] = number_geometrical_border;
> > }
>
> There are no add_facets method in MeshEditor. See comment above. However
> once that is created you can do something like:
>
> mesh.init(mesh.topology().dim()-1)
> FacetFunction facets(mesh);
> for(uint i = 0; i < num_side_edges; ++i)
> {
> Vertex v(mesh, a(i));
> for (FacetIterator(v) f; !f.end(), f++)
> {
> // Check vertex connectivity for facet (hard wired for edges...)
> // NOTE: Only need one check if we know that a(i) < b(i)
> if((f.entities(0)[0]==a(i) && f.entities(0)[1]==b(i)) ||\
> (f.entities(0)[0]==b(i) && f.entities(0)[1]==a(i)))
> {
> facet[f] = number_geometrical_border;
> break;
> }
> }
> }
>
Thanks a lot Johan !!
But another question arises: the FacetFunctions is an object related to
the mesh, so if I have stored the information about the label of the boundary
edges in a FacetFunctions and I do something like
refine(new_mesh, mesh)
on the new_mesh I have to initialize again all the labels ?
Is there any clever way in which I can keep the border label "alive"?
Marco V.
> > 2) When I read a function from a .xml file, is there any way in which I
> > can know
> > if something like boundary_markers, sub-domains, colours had been
> > defined on the
> > mesh??
>
> We do not store a MeshFunction in the mesh anymore. Instead we store
> MeshDomains with MeshValueCollections. These only store values for
> entities which has a boundary value. In a MeshValueCollection are the
> mesh entity stored using the cell number and the local mesh-entity
> number. This interface has also seen some recent changes in the
> development version and it has not fully stabilized.
>
> You can check if a mesh has MeshDomains by:
>
> mesh.domains().is_empty()
>
> Johan
>
> > Marco V.
> >
> >> An alternative, if you can live with calling Python from within the
> >> internals of your Octave package, would be to work with strings that
> >> are passed to Python. So you could do something like
> >>
> >> mesh = <some Octave data structure for a mesh like (p, e, t)>
> >> A = fem_matrix("dot(grad(u), grad(v))*dx, mesh, "P1")
> >>
> >> --
> >> Anders
> >> _______________________________________________
> >> fenics mailing list
> >> [email protected]
> >> http://fenicsproject.org/mailman/listinfo/fenics
> >
> >
> > _______________________________________________
> > fenics mailing list
> > [email protected]
> > http://fenicsproject.org/mailman/listinfo/fenics
> >
>
> _______________________________________________
> fenics mailing list
> [email protected]
> http://fenicsproject.org/mailman/listinfo/fenics
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
