Mesh refinement and the ability to transfer the data to the quadrature points of the new mesh on parallel::shared::triangulation.

87 views

Skip to first unread message

Alberto Salvadori

unread,

Apr 21, 2020, 12:03:39 PM4/21/20

to deal.II User Group

Dear community

I have been struggling in these days on the mesh refinement. I am encountering a problem that, so far, I just was unable to sort out.

Therefore, I wonder i f I can get some help.

Shortly: I want to refine my mesh for a vector problem in mechanics. After solving the problem on an initial grid, the solution is stored in the std::vector< double > this->accumulated_displacement

Moreover, I collected some further data in a constitutive class, MechanicalModels_FiniteDifference_Integrator<dim>* pmech;

associated to the cell->user_pointer()) pointer . On refinement, I aim at interpolating the solution and the data at the Gauss points, too.

This problem has been the subject of a few discussions and suggestions have been provided for parallel::distributed::triangulations. At present I am still using parallel::shared , though.

The step-26 shows very clearly how to deal with refinement and solution update. In fact, I copied the approach and it works very well.

The code gallery example 'Goal-oriented mesh adaptivity in elastoplasticity problems' seems to address the problem of how to propagate GP data. Being inspired by it, I wrote some code, that I am attaching below.

It turns out that it works well on 1 processor, but it fails in parallel. To test the code, I run a trivial test in which a uniform tensor

Fatq [ q_point ] = 1,0.15, 0, 0, 1.05, 0,0, 0.1,1

is stored at all GPs. After refinement on 1 processor, all GPs of the new triangulation have the same tensor Fatq. Once running on 4 processor, though, a print of Fatqs at processor 0 shows the following:

Problem LargeStrainMechanicsSolver_OneField_WithRemeshing defined

Reading material parameters from file ../input/mech_test_hex.materials ...

Reading refinement parameters from file ../input/mech_test_hex.msh_refinement ... done

Reading time discretization parameters from file ../input/mech_test_hex.time_discretization ... done

Time = 0.0000, step = 0

Initialization

Reading discretization from file ../input/mech_test_hex.msh ... done

Number of active cells: 8 (by partition: 2+2+2+2)

Number of degrees of freedom: 81 (by partition: 36+18+18+9)

Dirichlet faces: 24, Neumann faces (with non-zero tractions): 0, contact faces: 0

NR it. 0, Assembling..., convergence achieved.

Writing output..., 0.00 s. Elapsed time 0.02 s.

Time = 0.0500, step = 1

Refinement level: 0:

Number of active cells: 8 (by partition: 2+2+2+2)

Number of degrees of freedom: 81 (by partition: 36+18+18+9)

Dirichlet faces: 24, symmetry faces: 0

Dirichlet faces: 24, Neumann faces (with non-zero tractions): 0, contact faces: 0

NR it. 0, Assembling..., 0.00 s, norm of residual is 6034.853338302001248 Bicgstab , solver converged in 0 iterations, 0.01 s, updating q. p. data, 0.00 s.

NR it. 1, Assembling..., 0.00 s, norm of residual is 0.000000000000102 residual / initial_residual 0.000000000000000, convergence achieved.

Refinement level: 1:

Refining the grid, refined and coarsened fixed number, limited the refinement levels, executed coarsening and refinement, displacements transferred, quadrature_point_fields_trans interpolated,

Fatq [ q_point ] = 0, 0, 0, 0, 0, 0, 0, 0, 0

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.7347234759768e-18, 0, 1.05, -4.33680868994199e-19, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -2.31111593326468e-33, 0, 1.05, -3.05741378671474e-33, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 1, 0.15, -3.46944695195362e-18, 0, 1.05, -3.46944695195362e-18, 0, 0.1, 1

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.38777878078144e-17, 0, 1.05, -1.04083408558608e-17, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 1, 0.15, -8.67361737988401e-19, 0, 1.05, 1.10740971802266e-33, 0, 0.1, 1

Fatq [ q_point ] = 0.999999999999999, 0.15, 3.46944695195362e-18, 0, 1.05, 2.0222264416066e-33, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 1, 0.15, 3.46944695195362e-18, 0, 1.05, 1.73472347597681e-18, 0, 0.1, 1

Fatq [ q_point ] = 1, 0.15, -2.77555756156289e-17, 0, 1.05, -3.46944695195362e-18, 0, 0.1, 1

Fatq [ q_point ] = 0.999999999999999, 0.15, 4.69091588635142e-19, 0, 1.05, 1.52391206924327e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, 1.81732352981118e-18, 0, 1.05, 2.49138835915141e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -3.63711455892272e-18, 0, 1.05, -3.0838074433328e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -4.25316574819257e-18, 0, 1.05, -4.05868623114201e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, 2.6994902636017e-18, 0, 1.05, 1.02059236936956e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, 5.76972646078257e-18, 0, 1.05, 1.78043705753716e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -3.88766166070085e-18, 0, 1.05, -2.02424282585909e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -5.70905445967673e-18, 0, 1.05, -2.93570267228746e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -6.64306608528673e-18, 0, 1.05, -6.45689223356498e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -8.69707232649342e-18, 0, 1.05, -8.85367362460269e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.07492722328446e-17, 0, 1.05, -1.10646117461411e-17, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.47675616044972e-17, 0, 1.05, -1.54037482148961e-17, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -8.7097915464654e-18, 0, 1.05, -4.25321688944035e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.41121053024052e-17, 0, 1.05, -6.38815657011324e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 1, 0.15, -1.5296943470768e-17, 0, 1.05, -7.298052084669e-18, 0, 0.0999999999999999, 1

Fatq [ q_point ] = 0.999999999999999, 0.15, -2.55908862228645e-17, 0, 1.05, -1.11042962999379e-17, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, 4.33225541481152e-18, 0, 1.05, 6.52136776611694e-19, 0, 0.0999999999999998, 0.999999999999999

Fatq [ q_point ] = 0.999999999999998, 0.15, 8.66308621823777e-18, 0, 1.05, 1.25998458304831e-18, 0, 0.0999999999999999, 0.999999999999998

Fatq [ q_point ] = 0.999999999999999, 0.15, -4.07107486889156e-18, 0, 1.05, -1.24858769196604e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -6.77483896664912e-18, 0, 1.05, -2.11362165114141e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, 6.56265408977808e-18, 0, 1.05, 1.48817076737992e-19, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, 1.26154891492092e-17, 0, 1.05, 5.49033281434069e-19, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 1, 0.15, -4.32162197066969e-18, 0, 1.05, -1.89023074492321e-19, 0, 0.1, 1

Fatq [ q_point ] = 0.999999999999999, 0.15, -8.23072767813328e-18, 0, 1.05, -9.90638092286857e-19, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.02227395893444e-17, 0, 1.05, -2.64001457415429e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.80761845654335e-17, 0, 1.05, -4.5832728192994e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.86260698730475e-17, 0, 1.05, -4.54073904273203e-18, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 1, 0.15, -3.35141097503204e-17, 0, 1.05, -7.95687905348911e-18, 0, 0.1, 1

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.22894650505231e-17, 0, 1.05, -4.36339230029663e-19, 0, 0.0999999999999999, 0.999999999999999

Fatq [ q_point ] = 1, 0.15, -2.34912175413452e-17, 0, 1.05, -2.11775576480995e-18, 0, 0.1, 1

Fatq [ q_point ] = 1, 0.15, -2.31737411109709e-17, 0, 1.05, -7.74179381259976e-19, 0, 0.1, 1

Fatq [ q_point ] = 1, 0.15, -4.43374343686877e-17, 0, 1.05, -3.65742713853087e-18, 0, 0.1, 1

transferred the accumulated solution and quadrature points data, redifined the boundaries, done.

Number of active cells: 22 (by partition: 6+5+5+6)

Number of degrees of freedom: 180 (by partition: 78+33+42+27)

Dirichlet faces: 48, symmetry faces: 0

Dirichlet faces: 48, Neumann faces (with non-zero tractions): 0, contact faces: 0

NR it. 0, Assembling..., 0.01 s, norm of residual is 9210.109155179914524 Bicgstab , solver converged in 6 iterations, 0.00 s, updating q. p. data, 0.00 s.

NR it. 1, Assembling..., 0.01 s, norm of residual is 122.481867656167097 residual / initial_residual 0.013298633663563 Bicgstab , solver converged in 7 iterations, 0.00 s, updating q. p. data, 0.00 s.

NR it. 2, Assembling..., 0.01 s, norm of residual is 0.249221545837281 residual / initial_residual 0.000027059564837 Bicgstab , solver converged in 7 iterations, 0.00 s, updating q. p. data, 0.00 s.

NR it. 3, Assembling..., 0.01 s, norm of residual is 0.000001499773869 residual / initial_residual 0.000000000162840, convergence achieved.

Writing output..., 0.00 s. Elapsed time 0.06 s.

It seems therefore that one out of the 6 active cells on processor 0 is not providing the right answer, but I could not figure out the source of this issue. I surmise it is related to the parallel:shared triangulation, can any of you perhaps provide any suggestion?

I do appreciate your help very much

Kind regards

Alberto

template <int dim>

void LargeStrainMechanicsSolver_OneField_WithRemeshing<dim>

refine_the_grid ( bool neumann )

// this method refines the grid, by refining the triangulation and the

// quadrature point constitutive data.

// In a nutshell, parameters in the refine_and_coarsen_fixed_number

// are such that the overall amount of elements duplicates.

{

this->pcout << " Refining the grid, " << std::flush ;

// Make a field variable for history variables to be able to transfer the data to the quadrature points of the new mesh.

typedef ttl::Tensor<2, dim, double> FTensors;

FE_DGQ<dim> history_fe (1);

DoFHandler<dim> history_dof_handler ( this->triangulation );

history_dof_handler.distribute_dofs ( history_fe );

const unsigned int n_q_points = this->quadrature_formula.size();

std::vector< std::vector< Vector<double> > >

quadrature_point_fields_backup (dim, std::vector< Vector<double> >(dim)),

local_values_at_qpoints_backup (dim, std::vector< Vector<double> >(dim)),

local_fe_values_backup (dim, std::vector< Vector<double> >(dim));

for (unsigned i=0; i<dim; i++)

for (unsigned j=0; j<dim; j++)

{

quadrature_point_fields_backup[i][j].reinit( history_dof_handler.n_dofs() );

local_values_at_qpoints_backup[i][j].reinit( n_q_points );

local_fe_values_backup[i][j].reinit( history_fe.dofs_per_cell );

}

FullMatrix<double> qpoint_to_dof_matrix (

history_fe.dofs_per_cell,

n_q_points

);

FETools::compute_projection_from_quadrature_points_matrix(

history_fe,

this->quadrature_formula,

qpoint_to_dof_matrix

);

typename DoFHandler<dim>::active_cell_iterator

cell = this->dof_handler.begin_active() ,

endc = this->dof_handler.end() ,

dg_cell = history_dof_handler.begin_active();

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

if (cell->is_locally_owned())

{

typedef MechanicalModels_FiniteDifference_Integrator<dim>* pmech;

pmech *local_integrator = reinterpret_cast< pmech *>(cell->user_pointer());

Assert (local_integrator >=

&this->quadrature_point_integrators.front(),

ExcInternalError());

Assert (local_integrator <

&this->quadrature_point_integrators.back(),

ExcInternalError());

for (unsigned i=0; i<dim; i++)

for (unsigned j=0; j<dim; j++)

{

for (unsigned q_point=0; q_point<n_q_points; ++q_point)

{

FTensors Fatq = local_integrator[q_point]->get_Fnp1() ;

local_values_at_qpoints_backup[i][j]( q_point ) = Fatq( i,j );

}

qpoint_to_dof_matrix.vmult ( local_fe_values_backup[i][j] , local_values_at_qpoints_backup[i][j] );

dg_cell->set_dof_values ( local_fe_values_backup[i][j], quadrature_point_fields_backup[i][j] ) ;

}

// estimate local error via KellyErrorEstimator, refine and coarsen by fixed number

double fraction_of_cells = 0;

if ( dim == 2 ) fraction_of_cells = 1.0 / 3.0;

if ( dim == 3 ) fraction_of_cells = 1.0 / 7.0;

Vector<float> error_per_cell (this->triangulation.n_active_cells());

KellyErrorEstimator<dim>::estimate (

this->dof_handler,

QGauss<dim-1>( this->fe.degree + 1 ),

std::map<types::boundary_id, const Function<dim> *>(),

this->accumulated_displacement,

error_per_cell,

ComponentMask(),

MultithreadInfo::n_threads(),

this->this_mpi_process

);

const unsigned int n_local_cells = this->triangulation.n_locally_owned_active_cells ();

PETScWrappers::MPI::Vector

distributed_error_per_cell (this->mpi_communicator,

this->triangulation.n_active_cells(),

n_local_cells);

for (unsigned int i=0; i<error_per_cell.size(); ++i)

if (error_per_cell(i) != 0)

distributed_error_per_cell(i) = error_per_cell(i);

distributed_error_per_cell.compress (VectorOperation::insert);

error_per_cell = distributed_error_per_cell;

GridRefinement::refine_and_coarsen_fixed_number (this->triangulation,

error_per_cell,

fraction_of_cells, fraction_of_cells / 10.0 );

this->pcout << "refined and coarsened fixed number, " << std::flush ;

// limit the refinement levels

unsigned

min_grid_level = this->initial_global_refinements ,

max_grid_level = this->initial_global_refinements + this->n_adaptive_refinement_steps ;

if ( this->triangulation.n_levels() > max_grid_level)

{

for (const auto &cell :

this->triangulation.active_cell_iterators_on_level( max_grid_level ) )

cell->clear_refine_flag();

};

for (const auto &cell :

this->triangulation.active_cell_iterators_on_level( min_grid_level ))

cell->clear_coarsen_flag();

this->pcout << "limited the refinement levels, " << std::flush ;

this->triangulation.prepare_coarsening_and_refinement();

// define the data structures required for the

// transfer the accumulated solution and the constitutive data

// at Gauss points

this->accumulated_displacement_backup = this->accumulated_displacement ;

SolutionTransfer<dim> solution_trans( this->dof_handler );

solution_trans.prepare_for_coarsening_and_refinement( this->accumulated_displacement_backup ) ; // previous_solution );

SolutionTransfer< dim, Vector<double> > quadrature_point_fields_trans_0 ( history_dof_handler ) ;

quadrature_point_fields_trans_0.prepare_for_coarsening_and_refinement( quadrature_point_fields_backup[ 0 ] );

SolutionTransfer< dim, Vector<double> > quadrature_point_fields_trans_1 ( history_dof_handler ) ;

if (dim>1) quadrature_point_fields_trans_1.prepare_for_coarsening_and_refinement( quadrature_point_fields_backup[ 1 ] );

SolutionTransfer< dim, Vector<double> > quadrature_point_fields_trans_2 ( history_dof_handler ) ;

if (dim>2) quadrature_point_fields_trans_2.prepare_for_coarsening_and_refinement( quadrature_point_fields_backup[ 2 ] );

// execute coarsening and refinement

this->triangulation.execute_coarsening_and_refinement();

this->setup_system();

this->setup_quadrature_point_integrators ();

this->pcout << "executed coarsening and refinement, " << std::flush ;

// transfer the accumulated solution and the constitutive data

// at Gauss points via interpolation

this->accumulated_displacement.reinit ( this->dof_handler.n_dofs() );

solution_trans.interpolate( this->accumulated_displacement_backup, this->accumulated_displacement ); // previous_solution, solution);

this->hanging_node_constraints.distribute( this->accumulated_displacement ) ; // solution);

this->pcout << "displacements transferred, " << std::flush ;

history_dof_handler.distribute_dofs ( history_fe );

std::vector< std::vector< Vector<double> > >

quadrature_point_fields (dim, std::vector< Vector<double> >(dim)) ,

local_values_at_qpoints (dim, std::vector< Vector<double> >(dim)) ,

local_fe_values (dim, std::vector< Vector<double> >(dim)) ;

for (unsigned int i=0; i<dim; i++)

for (unsigned int j=0; j<dim; j++)

{

quadrature_point_fields[i][j].reinit( history_dof_handler.n_dofs() );

local_values_at_qpoints[i][j].reinit( n_q_points );

local_fe_values[i][j].reinit( history_fe.dofs_per_cell );

}

quadrature_point_fields_trans_0.interpolate( quadrature_point_fields_backup[ 0 ], quadrature_point_fields [ 0 ] );

if ( dim > 1) quadrature_point_fields_trans_1.interpolate( quadrature_point_fields_backup[ 1 ], quadrature_point_fields[ 1 ]);

if ( dim > 2) quadrature_point_fields_trans_2.interpolate( quadrature_point_fields_backup[ 2 ], quadrature_point_fields[ 2 ]);

this->pcout << "quadrature_point_fields_trans interpolated, " << std::flush ;

// Transfer the history data to the quadrature points of the new mesh

// In a final step, we have to get the data back from the now interpolated global

// field to the quadrature points on the new mesh. The following code will do that:

FullMatrix<double> dof_to_qpoint_matrix ( n_q_points, history_fe.dofs_per_cell );

FETools::compute_interpolation_to_quadrature_points_matrix(

history_fe,

this->quadrature_formula,

dof_to_qpoint_matrix

);

cell = this->dof_handler.begin_active();

endc = this->dof_handler.end();

dg_cell = history_dof_handler.begin_active();

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

if (cell->is_locally_owned())

{

typedef MechanicalModels_FiniteDifference_Integrator<dim>* pmech;

pmech *local_integrator = reinterpret_cast< pmech *>(cell->user_pointer());

Assert (local_integrator >=

&this->quadrature_point_integrators.front(),

ExcInternalError());

Assert (local_integrator <

&this->quadrature_point_integrators.back(),

ExcInternalError());

std::vector< FTensors > Fatq( n_q_points );

for (unsigned i=0; i<dim; i++)

for (unsigned j=0; j<dim; j++)

{

dg_cell->get_dof_values ( quadrature_point_fields[i][j], local_fe_values[i][j] );

dof_to_qpoint_matrix.vmult (local_values_at_qpoints[i][j], local_fe_values[i][j] );

for ( unsigned q_point=0; q_point<n_q_points; ++q_point )

Fatq[ q_point ]( i,j ) = local_values_at_qpoints[ i ][ j ]( q_point ) ;

}

// debug

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

for (unsigned q_point=0; q_point<n_q_points; ++q_point)

{

// debug

std::string Fstr = "";

PrintInVectorForm( Fatq [ q_point ], Fstr );

this->pcout << "Fatq [ q_point ] = " << Fstr << std::flush ;

if ( FrobeniusNorm( Fatq [ q_point ] ) > 0.1 ) // this should not be necessary

local_integrator[ q_point ]-> StepUpdate( Fatq [ q_point ], 0, false ) ;

}

};

this->pcout << " transferred the accumulated solution and quadrature points data, " << std::flush ;

// Re-Definition of boundaries

// symmcond was defined in create coarse grid and it is not updated here

// ----------------------------

this->pcout << " redifined the boundaries, done. \n" << std::flush ;

this->define_boundaries( this->symmcond, neumann );

}

Wolfgang Bangerth

unread,

Apr 21, 2020, 4:57:19 PM4/21/20

to dea...@googlegroups.com

Hi Alberto,
that's too much code for anyone to look through -- do you think you could come
up with a minimal, self-contained testcase that shows the problem? It's
entirely possible that there is a bug, but it would be very useful to have a
self-contained program that demonstrates it and that others can use for debugging.

Best
W.

On 4/21/20 10:03 AM, Alberto Salvadori wrote:
>
> Dear community
>
> I have been struggling in these days on the mesh refinement. I am encountering
> a problem that, so far, I just was unable to sort out.
> Therefore, I wonder i f I can get some help.
>
> Shortly: I want to refine my mesh for a vector problem in mechanics. After
> solving the problem on an initial grid, the solution is stored in the

> std::vector< double > *this*->accumulated_displacement

> *template*<*int*dim>
>
> *void*LargeStrainMechanicsSolver_OneField_WithRemeshing<dim>
>
> ::
>
> refine_the_grid( *bool*neumann )

>
>
> //
>
> // this method refines the grid, by refining the triangulation and the
>
> // quadrature point constitutive data.
>
> // In a nutshell, parameters in the refine_and_coarsen_fixed_number
>
> // are such that the overall amount of elements duplicates.
>
> //
>
>
> {
>

> *this*->pcout << " Refining the grid, " << std::flush ;

>
> //
>
> // Make a field variable for history variables to be able to transfer the data
> to the quadrature points of the new mesh.
>
> //
>

> *typedef* ttl::Tensor<2, dim, *double*> FTensors;
>
> FE_DGQ<dim> history_fe (1);
>
> DoFHandler<dim> history_dof_handler ( *this*->triangulation );
>
> history_dof_handler.distribute_dofs ( history_fe );
>
> *const* *unsigned* *int* n_q_points = *this*->quadrature_formula.size();
>
> std::vector< std::vector< Vector<*double*> > >
>
> quadrature_point_fields_backup (dim, std::vector< Vector<*double*> >(dim)),
>
> local_values_at_qpoints_backup (dim, std::vector< Vector<*double*> >(dim)),
>
> local_fe_values_backup (dim, std::vector< Vector<*double*> >(dim));
>
> *for* (*unsigned* i=0; i<dim; i++)
>
> *for* (*unsigned* j=0; j<dim; j++)

>
> {
>
> quadrature_point_fields_backup[i][j].reinit(
> history_dof_handler.n_dofs() );
>
> local_values_at_qpoints_backup[i][j].reinit( n_q_points );
>
> local_fe_values_backup[i][j].reinit( history_fe.dofs_per_cell );
>
> }
>

> FullMatrix<*double*> qpoint_to_dof_matrix (

>
>    history_fe.dofs_per_cell,
>
>    n_q_points
>
>    );
>
>
> FETools::compute_projection_from_quadrature_points_matrix(
>
> history_fe,
>

> *this*->quadrature_formula,
>
> *this*->quadrature_formula,
>
> qpoint_to_dof_matrix
>
> );
>
> *typename* DoFHandler<dim>::active_cell_iterator
>
> cell = *this*->dof_handler.begin_active() ,
>
> endc = *this*->dof_handler.end() ,
>
> dg_cell = history_dof_handler.begin_active();
>
> *for* (; cell!=endc; ++cell, ++dg_cell )
>
> *if* (cell->is_locally_owned())
>
> {
>
> *typedef*MechanicalModels_FiniteDifference_Integrator<dim>* pmech;
>
> pmech *local_integrator = *reinterpret_cast*< pmech *>(cell->user_pointer());
>
> Assert (local_integrator >=
>
> &*this*->quadrature_point_integrators.front(),
>
> ExcInternalError());
>
> Assert (local_integrator <
>
> &*this*->quadrature_point_integrators.back(),
>
> ExcInternalError());
>
> *for* (*unsigned* i=0; i<dim; i++)
>
> *for* (*unsigned* j=0; j<dim; j++)
>
> {
>
> *for* (*unsigned* q_point=0; q_point<n_q_points; ++q_point)

>
> {
>
> FTensors Fatq = local_integrator[q_point]->get_Fnp1() ;
>
> local_values_at_qpoints_backup[i][j]( q_point ) = Fatq( i,j );
>
> }
>
> qpoint_to_dof_matrix.vmult ( local_fe_values_backup[i][j] ,
> local_values_at_qpoints_backup[i][j] );
>
> dg_cell->set_dof_values ( local_fe_values_backup[i][j],
> quadrature_point_fields_backup[i][j] ) ;
>
> }
>
> }
>
> //
>
> // estimate local error via KellyErrorEstimator, refine and coarsen by fixed
> number
>
> //
>

> *double* fraction_of_cells = 0;
>
> *if* ( dim == 2 ) fraction_of_cells = 1.0 / 3.0;
>
> *if* ( dim == 3 ) fraction_of_cells = 1.0 / 7.0;
>
> Vector<*float*> error_per_cell (*this*->triangulation.n_active_cells());
>
> KellyErrorEstimator<dim>::estimate (
>
> *this*->dof_handler,
>
> QGauss<dim-1>( *this*->fe.degree + 1 ),
>
> std::map<types::boundary_id, *const* Function<dim> *>(),
>
> *this*->accumulated_displacement,

>
> error_per_cell,
>
> ComponentMask(),
>
> 0,
>
> MultithreadInfo::n_threads(),
>

> *this*->this_mpi_process
>
> );
>
> *const* *unsigned* *int* n_local_cells =
> *this*->triangulation.n_locally_owned_active_cells ();
>
> PETScWrappers::MPI::Vector
>
> distributed_error_per_cell (*this*->mpi_communicator,
>
> *this*->triangulation.n_active_cells(),
>
> n_local_cells);
>
> *for* (*unsigned* *int* i=0; i<error_per_cell.size(); ++i)
>
> *if* (error_per_cell(i) != 0)

>
> distributed_error_per_cell(i) = error_per_cell(i);
>
> distributed_error_per_cell.compress (VectorOperation::insert);
>
> error_per_cell = distributed_error_per_cell;
>

> GridRefinement::refine_and_coarsen_fixed_number (*this*->triangulation,

>
> error_per_cell,
>
> fraction_of_cells,
> fraction_of_cells / 10.0 );
>

> *this*->pcout << "refined and coarsened fixed number, "<< std::flush ;

>
> //
>
> // limit the refinement levels
>
> //
>

> *unsigned*
>
> min_grid_level = *this*->initial_global_refinements ,
>
> max_grid_level = *this*->initial_global_refinements +
> *this*->n_adaptive_refinement_steps ;
>
> *if* ( *this*->triangulation.n_levels() > max_grid_level)
>
> {
>
> *for* (*const* *auto* &cell :
>
> *this*->triangulation.active_cell_iterators_on_level( max_grid_level ) )
>
> cell->clear_refine_flag();
>
> };
>
> *for* (*const* *auto* &cell :
>
> *this*->triangulation.active_cell_iterators_on_level( min_grid_level ))
>
> cell->clear_coarsen_flag();
>
> *this*->pcout << "limited the refinement levels, "<< std::flush ;
>
> *this*->triangulation.prepare_coarsening_and_refinement();

>
> //
>
> // define the data structures required for the
>
> // transfer the accumulated solution and the constitutive data
>
> // at Gauss points
>
> //
>

> *this*->accumulated_displacement_backup = *this*->accumulated_displacement ;
>
> SolutionTransfer<dim> solution_trans( *this*->dof_handler );
>
> solution_trans.prepare_for_coarsening_and_refinement(
> *this*->accumulated_displacement_backup ) ; // previous_solution );
>
> SolutionTransfer< dim, Vector<*double*> > quadrature_point_fields_trans_0 (

> history_dof_handler ) ;
>
> quadrature_point_fields_trans_0.prepare_for_coarsening_and_refinement(
> quadrature_point_fields_backup[ 0 ] );
>

> SolutionTransfer< dim, Vector<*double*> > quadrature_point_fields_trans_1 (
> history_dof_handler ) ;
>
> *if* (dim>1)
> quadrature_point_fields_trans_1.prepare_for_coarsening_and_refinement(
> quadrature_point_fields_backup[ 1 ] );
>
> SolutionTransfer< dim, Vector<*double*> > quadrature_point_fields_trans_2 (
> history_dof_handler ) ;
>
> *if* (dim>2)

> quadrature_point_fields_trans_2.prepare_for_coarsening_and_refinement(
> quadrature_point_fields_backup[ 2 ] );
>
> //
>
> // execute coarsening and refinement
>
> //
>

> *this*->triangulation.execute_coarsening_and_refinement();
>
> *this*->setup_system();
>
> *this*->setup_quadrature_point_integrators ();
>
> *this*->pcout << "executed coarsening and refinement, "<< std::flush ;

>
> //
>
> // transfer the accumulated solution and the constitutive data
>
> // at Gauss points via interpolation
>
> //
>

> *this*->accumulated_displacement.reinit ( *this*->dof_handler.n_dofs() );
>
> solution_trans.interpolate( *this*->accumulated_displacement_backup,
> *this*->accumulated_displacement ); // previous_solution, solution);
>
> *this*->hanging_node_constraints.distribute( *this*->accumulated_displacement
> ) ; // solution);
>
> *this*->pcout << "displacements transferred, "<< std::flush ;
>
> history_dof_handler.distribute_dofs ( history_fe );
>
> std::vector< std::vector< Vector<*double*> > >
>
> quadrature_point_fields (dim, std::vector< Vector<*double*> >(dim)) ,
>
> local_values_at_qpoints (dim, std::vector< Vector<*double*> >(dim)) ,
>
> local_fe_values (dim, std::vector< Vector<*double*> >(dim)) ;
>
> *for* (*unsigned* *int* i=0; i<dim; i++)
>
> *for* (*unsigned* *int* j=0; j<dim; j++)

>
> {
>
> quadrature_point_fields[i][j].reinit( history_dof_handler.n_dofs() );
>
> local_values_at_qpoints[i][j].reinit( n_q_points );
>
> local_fe_values[i][j].reinit( history_fe.dofs_per_cell );
>
> }
>
> quadrature_point_fields_trans_0.interpolate(
> quadrature_point_fields_backup[ 0 ], quadrature_point_fields [ 0 ] );
>

> *if* ( dim > 1) quadrature_point_fields_trans_1.interpolate(

> quadrature_point_fields_backup[ 1 ], quadrature_point_fields[ 1 ]);
>

> *if* ( dim > 2) quadrature_point_fields_trans_2.interpolate(

> quadrature_point_fields_backup[ 2 ], quadrature_point_fields[ 2 ]);
>

> *this*->pcout << "quadrature_point_fields_trans interpolated, "<< std::flush ;

>
> //
>
> // Transfer the history data to the quadrature points of the new mesh
>
> // In a final step, we have to get the data back from the now interpolated global
>
> // field to the quadrature points on the new mesh. The following code will do
> that:
>
> //
>

> FullMatrix<*double*> dof_to_qpoint_matrix ( n_q_points,

> history_fe.dofs_per_cell );
>
> FETools::compute_interpolation_to_quadrature_points_matrix(
>
> history_fe,
>

> *this*->quadrature_formula,
>
> dof_to_qpoint_matrix
>
> );
>
> cell = *this*->dof_handler.begin_active();
>
> endc = *this*->dof_handler.end();
>
> dg_cell = history_dof_handler.begin_active();
>
> *for* (; cell != endc; ++cell, ++dg_cell)
>
> *if* (cell->is_locally_owned())
>
> {
>
> *typedef*MechanicalModels_FiniteDifference_Integrator<dim>* pmech;
>
> pmech *local_integrator = *reinterpret_cast*< pmech *>(cell->user_pointer());
>
> Assert (local_integrator >=
>
> &*this*->quadrature_point_integrators.front(),
>
> ExcInternalError());
>
> Assert (local_integrator <
>
> &*this*->quadrature_point_integrators.back(),

>
> ExcInternalError());
>
> std::vector< FTensors > Fatq( n_q_points );
>

> *for* (*unsigned* i=0; i<dim; i++)
>
> *for* (*unsigned* j=0; j<dim; j++)

>
> {
>
> dg_cell->get_dof_values ( quadrature_point_fields[i][j], local_fe_values[i][j] );
>
> dof_to_qpoint_matrix.vmult (local_values_at_qpoints[i][j],
> local_fe_values[i][j] );
>

> *for* ( *unsigned* q_point=0; q_point<n_q_points; ++q_point )

>
> Fatq[ q_point ]( i,j ) = local_values_at_qpoints[ i ][ j ](
> q_point ) ;
>
> }
>
> // debug
>

> *this*->pcout << "\n" << std::flush ;
>
> *for* (*unsigned* q_point=0; q_point<n_q_points; ++q_point)

>
> {
>
> // debug
>
> std::string Fstr = "";
>
> PrintInVectorForm( Fatq [ q_point ], Fstr );
>

> *this*->pcout << "Fatq [ q_point ] = " << Fstr << std::flush ;
>
> *if* ( FrobeniusNorm( Fatq [ q_point ] ) > 0.1 ) // this should not be necessary

>
> local_integrator[ q_point ]-> StepUpdate( Fatq [ q_point ], 0,

> *false* ) ;
>
> }
>
> };
>
> *this*->pcout << " transferred the accumulated solution and quadrature points

> data, "<< std::flush ;
>
> // Re-Definition of boundaries
>
> // symmcond was defined in create coarse grid and it is not updated here
>
> // ----------------------------
>

> *this*->pcout << " redifined the boundaries, done. \n"<< std::flush ;
>
> *this*->define_boundaries( *this*->symmcond, neumann );
>
> }
>
>
>
> Informativa sulla Privacy: http://www.unibs.it/node/8155
> <https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.unibs.it%2Fnode%2F8155&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7C8f6ac5c549b24d323bd508d7e60d914f%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637230818272783700&sdata=q7xhP6h53c57FVagGYS0oXfLph8rFvw42ST4tgIq514%3D&reserved=0>
>
> --
> The deal.II project is located at http://www.dealii.org/
> <https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7C8f6ac5c549b24d323bd508d7e60d914f%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637230818272793705&sdata=yN3FifMdpAXFEJBj%2BA6k%2FT1jNoSIIR%2BXtHrRWpBebzE%3D&reserved=0>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7C8f6ac5c549b24d323bd508d7e60d914f%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637230818272803690&sdata=2KHUeYqFnxd%2FXLz7Pv6y3FzVAlIQWYZ38Wod4dlRfcU%3D&reserved=0>
> ---
> 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
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/7629f35a-50e3-49f3-8fcd-c5370df83937%40googlegroups.com
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2F7629f35a-50e3-49f3-8fcd-c5370df83937%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=02%7C01%7CWolfgang.Bangerth%40colostate.edu%7C8f6ac5c549b24d323bd508d7e60d914f%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637230818272803690&sdata=eopX1HsS6HGH98x3GGbS%2FU8yS%2F00cFmDqATyH9lWzAg%3D&reserved=0>.

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

Alberto Salvadori

unread,

Apr 22, 2020, 4:36:15 AM4/22/20

to deal.II User Group

Sure, Wolfgang. I will do my best below.

To test the problem, one could move from step 18 or 26 and the method "refine_initial_grid" therein. Although some local data should be stored in the cell->user_pointer() , I am storing the deformation gradient tensor, to test the code below it is not necessary.

My test begins creating a cube

GridGenerator::hyper_cube (this->triangulation, -1, 1);

this->triangulation.refine_global ( 1 );

and the mesh looks consistent on 4 processors:

Number of active cells: 8 (by partition: 2+2+2+2)

Number of degrees of freedom: 81 (by partition: 36+18+18+9)

A trivial set of boundary conditions and body forces allow seeking for a manufactured solution, with deformation gradient

F = 1,0.15, 0, 0, 1.05, 0,0, 0.1,1

independent upon the position (i.e. uniform). This is likely not a good test for refinement, at any rate the method “refine_initial_grid" is invoked and indeed the mesh is refined:

Number of active cells: 22 (by partition: 6+5+5+6)

Number of degrees of freedom: 180 (by partition: 78+33+42+27)

The transfer of the nodal values solution proceeds smoothly following the approach in step-26, so I attacked the transfer of data stored at each GP (the deformation gradient) by the cell->user_pointer(). The code gallery example 'Goal-oriented mesh adaptivity in elastoplasticity problems' seems to address the problem of how to propagate GP data and corresponds to a few discussions on this group. Being inspired by it, I wrote some code, that follows.

// 1 - read data from the cell->user_pointer() and store them temporarily into

// the container quadrature_point_fields_backup

FE_DGQ<dim> history_fe (1);

DoFHandler<dim> history_dof_handler ( this->triangulation );

history_dof_handler.distribute_dofs ( history_fe );

const unsigned int n_q_points = this->quadrature_formula.size();

std::vector< std::vector< Vector<double> > >

quadrature_point_fields_backup (dim, std::vector< Vector<double> >(dim)),

local_values_at_qpoints_backup (dim, std::vector< Vector<double> >(dim)),

local_fe_values_backup (dim, std::vector< Vector<double> >(dim));

for (unsigned i=0; i<dim; i++)

for (unsigned j=0; j<dim; j++)

{

quadrature_point_fields_backup[i][j].reinit( history_dof_handler.n_dofs() );

local_values_at_qpoints_backup[i][j].reinit( n_q_points );

local_fe_values_backup[i][j].reinit( history_fe.dofs_per_cell );

}

FullMatrix<double> qpoint_to_dof_matrix (

history_fe.dofs_per_cell,

n_q_points

);

qpoint_to_dof_matrix *= 0;

FETools::compute_projection_from_quadrature_points_matrix(

history_fe,

this->quadrature_formula,

qpoint_to_dof_matrix

);

typename DoFHandler<dim>::active_cell_iterator

cell = this->dof_handler.begin_active() ,

endc = this->dof_handler.end() ,

dg_cell = history_dof_handler.begin_active();

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

if (cell->is_locally_owned())

{

for (unsigned i=0; i<dim; i++)

for (unsigned j=0; j<dim; j++)

{

for (unsigned q_point=0; q_point<n_q_points; ++q_point)

{

... read data at a GP from the cell->user_pointer() and store them in a tensor, Fatq ...

// in fact, to run a test case it is not really necessary to read Fatq, one could just define it

local_values_at_qpoints_backup[i][j]( q_point ) = Fatq( i,j );

}

qpoint_to_dof_matrix.vmult ( local_fe_values_backup[i][j] , local_values_at_qpoints_backup[i][j] );

dg_cell->set_dof_values ( local_fe_values_backup[i][j], quadrature_point_fields_backup[i][j] ) ;

}

The vector quadrature_point_fields_backup should at this point contain a copy of all F tensors at GPs for the locally owned cells, and indeed it happens. Estimation of error via KellyErrorEstimator is exactly as in step-26, thus I built the solution transfer and proceeded further:

// 3 - define the data structures required for the SolutionTransfer

SolutionTransfer< dim, Vector<double> > quadrature_point_fields_trans_0 ( history_dof_handler ) ;

quadrature_point_fields_trans_0.prepare_for_coarsening_and_refinement( quadrature_point_fields_backup[ 0 ] );

SolutionTransfer< dim, Vector<double> > quadrature_point_fields_trans_1 ( history_dof_handler ) ;

if (dim>1) quadrature_point_fields_trans_1.prepare_for_coarsening_and_refinement( quadrature_point_fields_backup[ 1 ] );

SolutionTransfer< dim, Vector<double> > quadrature_point_fields_trans_2 ( history_dof_handler ) ;

if (dim>2) quadrature_point_fields_trans_2.prepare_for_coarsening_and_refinement( quadrature_point_fields_backup[ 2 ] );

// execute coarsening and refinement

this->triangulation.execute_coarsening_and_refinement();

this->setup_system();

this->setup_quadrature_point_integrators ();

this->pcout << "executed coarsening and refinement, " << std::flush ;

// 4 - transfer the data back into the cell->user_pointer() data structure

// transfer the accumulated solution and the constitutive data

// at Gauss points via interpolation

history_dof_handler.distribute_dofs ( history_fe );

std::vector< std::vector< Vector<double> > >

quadrature_point_fields (dim, std::vector< Vector<double> >(dim)) ,

local_values_at_qpoints (dim, std::vector< Vector<double> >(dim)) ,

local_fe_values (dim, std::vector< Vector<double> >(dim)) ;

for (unsigned int i=0; i<dim; i++)

for (unsigned int j=0; j<dim; j++)

{

quadrature_point_fields[i][j].reinit( history_dof_handler.n_dofs() );

local_values_at_qpoints[i][j].reinit( n_q_points );

local_fe_values[i][j].reinit( history_fe.dofs_per_cell );

}

quadrature_point_fields_trans_0.interpolate( quadrature_point_fields_backup[ 0 ], quadrature_point_fields [ 0 ] );

if ( dim > 1) quadrature_point_fields_trans_1.interpolate( quadrature_point_fields_backup[ 1 ], quadrature_point_fields[ 1 ]);

if ( dim > 2) quadrature_point_fields_trans_2.interpolate( quadrature_point_fields_backup[ 2 ], quadrature_point_fields[ 2 ]);

// Transfer the history data to the quadrature points of the new mesh

// In a final step, we have to get the data back from the now interpolated global

// field to the quadrature points on the new mesh. The following code will do that:

FullMatrix<double> dof_to_qpoint_matrix ( n_q_points, history_fe.dofs_per_cell );

FETools::compute_interpolation_to_quadrature_points_matrix(

history_fe,

this->quadrature_formula,

dof_to_qpoint_matrix

);

cell = this->dof_handler.begin_active();

endc = this->dof_handler.end();

dg_cell = history_dof_handler.begin_active();

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

if (cell->is_locally_owned())

{

... define the std::vector of tensors Fatq( n_q_points ) ...

for (unsigned i=0; i<dim; i++)

for (unsigned j=0; j<dim; j++)

{

dg_cell->get_dof_values ( quadrature_point_fields[i][j], local_fe_values[i][j] );

dof_to_qpoint_matrix.vmult (local_values_at_qpoints[i][j], local_fe_values[i][j] );

for ( unsigned q_point=0; q_point<n_q_points; ++q_point )

Fatq[ q_point ]( i,j ) = local_values_at_qpoints[ i ][ j ]( q_point ) ;

}

for (unsigned q_point=0; q_point<n_q_points; ++q_point)

{

... write Fatq back into the cell->user_pointer() ...

}

};

Since the original F tensor is uniform, I would expect that after interpolation the set of F at GPs would be uniform as well. In fact, is it like that when running on 1 single processor. Once I attempted to run on 4 processors, at the node 0 there are 6 cells

Number of active cells: 22 (by partition: 6+5+5+6)

and the outcome of F is the following at a generic GP at each of the 6 cells:

cell 0:

Fatq [ q_point ] = 0, 0, 0, 0, 0, 0, 0, 0, 0

at all other 5 cells:

Fatq [ q_point ] = 0.999999999999999, 0.15, -1.7347234759768e-18, 0, 1.05, -4.33680868994199e-19, 0, 0.0999999999999999, 0.999999999999999

In conclusion, therefore, something goes wrong with the first cell. I assume I coded in an improper way but so far I was unable to find out the bug.

It turned out that filling completely the vector quadrature_point_fields_backup as follows:

double tmp = Utilities::MPI::sum( quadrature_point_fields_backup[i][j][q_point], this->mpi_communicator) ;

quadrature_point_fields_backup[i][j][q_point] = tmp;

sorts the problem out, but I am quite sure it should not be necessary.

Thanks for your help.

Alberto

> email to dea...@googlegroups.com
> <mailto:dea...@googlegroups.com>.

Marc Fehling

unread,

Apr 22, 2020, 12:17:32 PM4/22/20

to deal.II User Group

Hi Alberto!

If I understood you correctly, you transfer quadrature point data, with the `SolutionTransfer` class which is meant to transfer finite element approximations.

A different class dedicated to the transfer of quadrature point data already exists: It is called `TransferableQuadraturePointData`. Examples on how to use that feature can be found in `tests/base/quadrature_point_data_{02|03|04}`.

You could also use the `CellDataTransfer` class to transfer cell related data, i.e. stored as `Vector<Vector<double>>` in your case if I interpreted your code correctly. However, this particular feature is only available in the current `master` branch of the library and has not been released yet.

Hope this gives you some more options to find a solution!

Best,
Marc

Reply all

Reply to author

Forward

0 new messages