Extracting volume data to the boundary

Alberto Salvadori

unread,

Sep 3, 2019, 12:36:59 PM9/3/19

to deal.II User Group

Dear community

I am trying to find a good way to deal with the following problem.

I am interested in a Laplace-Beltrami simulation (say with unknown scalar field c) on a surface that advects (with displacement field u) enclosesing a volume of say elastic material.

The two problems are per se uncoupled, but for the geometrical evolution. In this regard, one may solve separately for displacements in the volume and for c on the surface, once the latter is informed on the displacement field (and related gradients) on the surface nodes (and Gauss points).

Extracting the surface mesh from the volume is quite simple. Which is though the best way to extract (project) the displacement solution on the boundary? Any suggestion/ hint?

Many thanks,

Alberto

Alberto Salvadori

unread,

Sep 7, 2019, 6:09:30 AM9/7/19

to deal.II User Group

Dear community

I apologize if the question was too generic. I attempt here to be more precise, asking for some suggestions.

In solving a Laplace-Beltrami problem on an advecting surface one should identify the Gauss points on the manifold dim-1 cell and retrieve at such locations relevant information from the solution of the advection problem, using the dim dof_handler of the volume. I wonder how to connect a manifold dim-1 cell to the volumetric cell it was extracted from. The GridGenerator::extract_boundary_mesh method seem to provide some information, since "it returns A map that for each cell of the surface mesh (key) returns an iterator to the corresponding face of a cell of the volume mesh (value). " . I am not sure how I can use this map to link a manifold cell to the corresponding volume cell.

Any help is very much appreciated. Thanks!

Alberto

Alberto Salvadori

unread,

Sep 9, 2019, 6:08:15 AM9/9/19

to deal.II User Group

Dear community

in case of use, I attempted to write a small code that seems to solve the issue of linking manifold cells to volume cells.

As first, two maps are defined:

// maps

// the map surface_to_volumefaces_mapping contains the outcome of the method extract_boundary_mesh and connects

// panels from the manifold triangulation to the faces on the boundary of the volume triangulation.

// the map manifoldcells_to_volumecells_mapping pairs surface panels from the manifold triangulation

// to the cells of the volume triangulation. It is used in the method pair_the_triangulations().

//

std::map<

typename parallel::shared::Triangulation<manifold_dim,dim>::cell_iterator,

typename parallel::shared::Triangulation<dim>::face_iterator

>

surface_to_volumefaces_mapping;

std::map<

typename parallel::shared::Triangulation<manifold_dim,dim>::cell_iterator,

typename parallel::shared::Triangulation<dim>::cell_iterator

>

manifoldcells_to_volumecells_mapping;

then the manifold grid is generated

//

// Surface triangulation definition

//

this->pcout << "\n Generation of surface discretization ... " << std::flush;

surface_to_volumefaces_mapping = GridGenerator::extract_boundary_mesh ( this->triangulation, this->manifold_triangulation );

this->manifold_triangulation.set_all_manifold_ids(0);

this->manifold_dof_handler.distribute_dofs ( this->manifold_fe );

this->manifold_accumulated_solution.reinit ( this->manifold_dof_handler.n_dofs() );

this->manifold_former_step_solution.reinit ( this->manifold_dof_handler.n_dofs() );

this->manifold_initial_guess.reinit ( this->manifold_dof_handler.n_dofs() );

this->pcout << " done \n" << std::flush;

//

// Surface triangulation and volume triangulation pairing

//

pair_the_triangulations();

In the method pair_the_triangulations() a manifold cell is linked to the corresponding volume cell for parallel::shared::triangulations as follows.

pair_the_triangulations ( bool write_map_report )

//

// this function associates each

// surface panel to the volume cell from which it emanates.

// the map manifoldcells_to_volumecells_mapping contains the pair surface/volume panels

//

{

// unsigned face_counter = 0;

this->pcout << "pairing the two triangulations ... " << std::flush;

// loop through the volumetric triangulation

typename parallel::shared::Triangulation<dim>::active_cell_iterator

cell = this->triangulation.begin_active (),

endc = this->triangulation.end();

for (; cell!=endc; ++cell)

// We do not seek for locally_owned on the volumetric triangulation

{

// debug

// this->pcout << " new cell " << cell->id() << " " << std::flush;

// select a face at boundary

for (unsigned int face_number=0; face_number<GeometryInfo<dim>::faces_per_cell; ++face_number)

{

const typename parallel::shared::Triangulation<dim>::face_iterator

face = cell->face( face_number );

if ( face->at_boundary() )

{

// scan the surface_to_volumefaces_mapping and

// associate maifold panels to volume panels

typename

std::map< typename parallel::shared::Triangulation<manifold_dim,dim>::cell_iterator,

typename parallel::shared::Triangulation<dim>::face_iterator >

::iterator

map_it = surface_to_volumefaces_mapping.begin(),

map_end = surface_to_volumefaces_mapping.end() ;

for (; map_it!=map_end; ++map_it)

if ( face == map_it->second )

{

// debug

// this->pcout << cell->face_index(face_number) << ", " << std::flush;

// face_counter++;

manifoldcells_to_volumecells_mapping[ map_it->first ] = cell ;

}

// debug

// this->pcout << "\n I found " << face_counter << " faces. \n bye. \n " << std::flush;

// map report

if ( write_map_report )

{

typename std::map< typename parallel::shared::Triangulation<manifold_dim,dim>::cell_iterator,

typename parallel::shared::Triangulation<dim>::cell_iterator >

::iterator

manifoldcells_to_volumecells_it = manifoldcells_to_volumecells_mapping.begin() ,

manifoldcells_to_volumecells_end = manifoldcells_to_volumecells_mapping.end() ;

this->pcout << "\n map report: " << std::flush;

for (; manifoldcells_to_volumecells_it!=manifoldcells_to_volumecells_end; ++manifoldcells_to_volumecells_it)

this->pcout << "[ " << manifoldcells_to_volumecells_it->first->id() << ", " << manifoldcells_to_volumecells_it->second->id() << " ], " << std::flush;

this->pcout << "\n done. \n" << std::flush;

}

// debug

// std::exit( 0 );

}

Comments and suggestions are very welcome.

Thank you in advance

Alberto

Wolfgang Bangerth

unread,

Sep 9, 2019, 1:09:09 PM9/9/19

to dea...@googlegroups.com

> In solving a Laplace-Beltrami problem on an advecting surface one should
> identify the Gauss points on the manifold dim-1 cell and retrieve at
> such locations relevant information from the solution of the advection
> problem, using the dim dof_handler of the volume. I wonder how to
> connect a manifold dim-1 cell to the volumetric cell it was extracted
> from. The GridGenerator::extract_boundary_mesh method seem to provide
> some information, since "it returns A map that for each cell of the
> surface mesh (key) returns an iterator to the corresponding face of a
> cell of the volume mesh (value). " . I am not sure how I can use this
> map to link a manifold cell to the corresponding volume cell.

Alberto,
the way this is supposed to be used is as follows:

When you are assembling the linear system for the surface problem, you
will have an FEValues<dim-1,dim> object for the surface shape functions.
You initialize it with a Quadrature<dim-1>.

But then you also need to evaluate the volume solution on the same
quadrature points. You do this by creating an FEFaceValues<dim> object,
which requires you to also provide a Quadrature<dim-1> object -- which
you want to choose the same as above.

Now, if you need the volume solution when assembling something for the
surface problem, you are sitting on a cell of the surface mesh; use the
map returned by the extract_boundary_mesh() function to obtain what the
corresponding face of the volume mesh is and reinit() the FEFaceValues
for that cell. You can then use the FEFaceValues object to evaluate the
volume solution at the same quadrature points as you use for the
FEValues object of the surface mesh. The only thing you may have to pay
attention to is that the surface cell may be inverted compared to the
volume cell's face -- in which case the quadrature points of the two
objects match, but are in a different order. You'll have to translate
between these permutations then.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

Alberto Salvadori

unread,

Sep 10, 2019, 6:15:49 AM9/10/19

to dea...@googlegroups.com

Thank you very much, Wolfgang, crystal clear as always.

It is still unclear to me, though, how to reinit() the FEFaceValues
for the cell making use of the

map returned by the extract_boundary_mesh() function. I learned to reinit

FEFaceValues as:

fe_face_values.reinit (cell, face_number);

Is there a "simple way" to get the (volume) cell iterator and the face_number from the

parallel::shared::Triangulation<dim>::face_iterator in the map returned by the extract_boundary_mesh() function?

Or shall I use a different way to reinit FEFaceValues?

Thank you

Alberto

Alberto Salvadori
Dipartimento di Ingegneria Meccanica e Industriale (DIMI)
Universita` di Brescia, via Branze 43, 25123 Brescia
Italy
tel 030 3711239
fax 030 3711312

e-mail:
alberto....@unibs.it
web-page:
http://m4lab.unibs.it/faculty.html

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a60b9b5e-fef3-fd88-5e54-65bed698a5bc%40colostate.edu.

Informativa sulla Privacy: http://www.unibs.it/node/8155

luca.heltai

unread,

Sep 11, 2019, 2:38:45 AM9/11/19

to Deal.II Users

Alberto,

be aware of the fact that, when you extract the boundary mesh, the orientation of the cells may be different w.r.t. to the corresponding face on the volumetric mesh. In other words, using

FEValues<dim, spacedim> codim_one(…);
FEFaceValues<dim> codim_zero(…);

codim_one.reinit(codim_one_cell);
codim_zero.reinit(codmi_zero_cell, corresponding_face);

in general, you will *NOT* get the same ordering of the quadrature points, and of the dof indices.

In the past, we proceeded by constructing a dof_map between codim zero and codim one (instead of trying to assemble together fe values coming from grids with different dimensions), or by constructing (by hand) the permutation of the dof indices.

As far as your question is concerned:

c1_to_c0 = GridGenerator::extract_boundary_mesh(triangulation, c1_tria);

// map volume mesh face -> codimension 1 mesh cell
for (auto c1_cell : c1_tria.active_cell_iterators())
c0_to_c1[c1_to_c0[c1_cell]] = c1_cell;

// generate a mapping that maps codimension-1 cells
// to codimension-0 cells and faces
for (auto cell :
triangulation
.active_cell_iterators()) // disp_dof.active_cell_iterators())
for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->face(f)->at_boundary())
{
auto c1_cell = c0_to_c1[cell->face(f)];
c1_to_c0_cells_and_faces[c1_cell] = {cell, f};
}

L.

Alberto Salvadori

unread,

Sep 15, 2019, 2:18:23 AM9/15/19

to dea...@googlegroups.com

Thank you very much, Luca.

Alberto Salvadori
Associate Professor
DIMI, University of Brescia, Italy

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/1DD9B280-CD2F-4C83-81FE-3FC00A522B91%40gmail.com.

--