Conversion of the LDG-example from the code gallery for using MeshWorker
20 views
Skip to first unread message
Maxi Miller
Mar 4, 2019, 8:53:45 AM3/4/19
to deal.II User Group
I tried to implement the LDG-example from the code gallery using MeshWorker::mesh_loop to get a better feeling for this function. For that I created three new functions, local_assemble_system(), local_assemble_boundaries() and local_assemble_faces(), and a fourth function for copying everything together. For the functions I simply copied the existing code (and modified it, when necessary):
template <int dim>
void LDGPoissonProblem<dim>::local_assemble_system(const typename DoFHandler<dim>::active_cell_iterator &cell,
AssemblyScratchData &scratch,
AssemblyCopyData &data)
{
scratch.fe_values.reinit(cell);
const unsigned int dofs_per_cell = scratch.fe_values.get_fe().dofs_per_cell;
const unsigned int n_q_points = scratch.fe_values.n_quadrature_points;
const FEValuesExtractors::Vector VectorField(0);
const FEValuesExtractors::Scalar Potential(dim);
std::vector<double> rhs_values(n_q_points);
rhs_function.value_list(scratch.fe_values.get_quadrature_points(),
rhs_values);
data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.vi_ui_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.vi_ue_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.ve_ui_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.ve_ue_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.cell_rhs.reinit(dofs_per_cell);
data.local_dof_indices.resize(dofs_per_cell);
data.cell_matrix = 0;
data.cell_rhs = 0;
data.vi_ui_matrix = 0;
data.vi_ue_matrix = 0;
data.ve_ui_matrix = 0;
data.ve_ue_matrix = 0;
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const Tensor<1, dim> psi_i_field = scratch.fe_values[VectorField].value(i,q);
const double div_psi_i_field = scratch.fe_values[VectorField].divergence(i,q);
const double psi_i_potential = scratch.fe_values[Potential].value(i,q);
const Tensor<1, dim> grad_psi_i_potential = scratch.fe_values[Potential].gradient(i,q);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const Tensor<1, dim> psi_j_field = scratch.fe_values[VectorField].value(j,q);
const double psi_j_potential = scratch.fe_values[Potential].value(j,q);
data.cell_matrix(i,j) += ( (psi_i_field * psi_j_field)
-
(div_psi_i_field * psi_j_potential)
-
(grad_psi_i_potential * psi_j_field)
) * scratch.fe_values.JxW(q);
}
data.cell_rhs(i) += psi_i_potential *
rhs_values[q] *
scratch.fe_values.JxW(q);
}
}
cell->get_dof_indices(data.local_dof_indices);
}
template <int dim>
void LDGPoissonProblem<dim>::local_assemble_faces(const typename DoFHandler<dim>::active_cell_iterator &interior_cell,
const unsigned int interior_face_index,
const unsigned int interior_subface_index,
const typename DoFHandler<dim>::active_cell_iterator &exterior_cell,
const unsigned int exterior_face_index,
const unsigned int exterior_subface_index,
AssemblyScratchData &scratch,
AssemblyCopyData &data)
{
if(interior_subface_index == numbers::invalid_unsigned_int)
{
double h = std::min(interior_cell->diameter(),
exterior_cell->diameter());
scratch.fe_face_values.reinit(interior_cell, interior_face_index);
scratch.fe_neighbour_face_values.reinit(exterior_cell, exterior_face_index);
const unsigned int n_face_points = scratch.fe_face_values.n_quadrature_points;
const unsigned int dofs_this_cell = scratch.fe_face_values.dofs_per_cell;
const unsigned int dofs_neighbor_cell = scratch.fe_neighbour_face_values.dofs_per_cell;
data.cell_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.cell_rhs.reinit(dofs_this_cell);
data.local_dof_indices.resize(dofs_this_cell);
data.local_neighbor_dof_indices.resize(dofs_neighbor_cell);
interior_cell->get_dof_indices(data.local_dof_indices);
exterior_cell->get_dof_indices(data.local_neighbor_dof_indices);
data.vi_ui_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.vi_ue_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.ve_ui_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.ve_ue_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.cell_matrix = 0;
data.cell_rhs = 0;
data.vi_ui_matrix = 0;
data.vi_ue_matrix = 0;
data.ve_ui_matrix = 0;
data.ve_ue_matrix = 0;
const FEValuesExtractors::Vector VectorField(0);
const FEValuesExtractors::Scalar Potential(dim);
Point<dim> beta;
for (int i=0; i<dim; i++)
beta(i) = 1.0;
beta /= sqrt(beta.square() );
for (unsigned int q=0; q<n_face_points; q++)
{
for (unsigned int i=0; i<dofs_this_cell; i++)
{
const Tensor<1,dim> psi_i_field_minus =
scratch.fe_face_values[VectorField].value(i,q);
const double psi_i_potential_minus =
scratch.fe_face_values[Potential].value(i,q);
for (unsigned int j=0; j<dofs_this_cell; j++)
{
const Tensor<1,dim> psi_j_field_minus =
scratch.fe_face_values[VectorField].value(j,q);
const double psi_j_potential_minus =
scratch.fe_face_values[Potential].value(j,q);
data.vi_ui_matrix(i,j) += (0.5 * (
psi_i_field_minus *
scratch.fe_face_values.normal_vector(q) *
psi_j_potential_minus
+
psi_i_potential_minus *
scratch.fe_face_values.normal_vector(q) *
psi_j_field_minus )
+
beta *
psi_i_field_minus *
psi_j_potential_minus
-
beta *
psi_i_potential_minus *
psi_j_field_minus
+
(scratch.penalty/h) *
psi_i_potential_minus *
psi_j_potential_minus
) *
scratch.fe_face_values.JxW(q);
}
for (unsigned int j=0; j<dofs_neighbor_cell; j++)
{
const Tensor<1,dim> psi_j_field_plus =
scratch.fe_neighbour_face_values[VectorField].value(j,q);
const double psi_j_potential_plus =
scratch.fe_neighbour_face_values[Potential].value(j,q);
data.vi_ue_matrix(i,j) += ( 0.5 * (
psi_i_field_minus *
scratch.fe_face_values.normal_vector(q) *
psi_j_potential_plus
+
psi_i_potential_minus *
scratch.fe_face_values.normal_vector(q) *
psi_j_field_plus )
-
beta *
psi_i_field_minus *
psi_j_potential_plus
+
beta *
psi_i_potential_minus *
psi_j_field_plus
-
(scratch.penalty/h) *
psi_i_potential_minus *
psi_j_potential_plus
) *
scratch.fe_face_values.JxW(q);
}
}
for (unsigned int i=0; i<dofs_neighbor_cell; i++)
{
const Tensor<1,dim> psi_i_field_plus =
scratch.fe_neighbour_face_values[VectorField].value(i,q);
const double psi_i_potential_plus =
scratch.fe_neighbour_face_values[Potential].value(i,q);
for (unsigned int j=0; j<dofs_this_cell; j++)
{
const Tensor<1,dim> psi_j_field_minus =
scratch.fe_face_values[VectorField].value(j,q);
const double psi_j_potential_minus =
scratch.fe_face_values[Potential].value(j,q);
data.ve_ui_matrix(i,j) += (-0.5 * (
psi_i_field_plus *
scratch.fe_face_values.normal_vector(q) *
psi_j_potential_minus
+
psi_i_potential_plus *
scratch.fe_face_values.normal_vector(q) *
psi_j_field_minus)
-
beta *
psi_i_field_plus *
psi_j_potential_minus
+
beta *
psi_i_potential_plus *
psi_j_field_minus
-
(scratch.penalty/h) *
psi_i_potential_plus *
psi_j_potential_minus
) *
scratch.fe_face_values.JxW(q);
}
for (unsigned int j=0; j<dofs_neighbor_cell; j++)
{
const Tensor<1,dim> psi_j_field_plus =
scratch.fe_neighbour_face_values[VectorField].value(j,q);
const double psi_j_potential_plus =
scratch.fe_neighbour_face_values[Potential].value(j,q);
data.ve_ue_matrix(i,j) += (-0.5 * (
psi_i_field_plus *
scratch.fe_face_values.normal_vector(q) *
psi_j_potential_plus
+
psi_i_potential_plus *
scratch.fe_face_values.normal_vector(q) *
psi_j_field_plus )
+
beta *
psi_i_field_plus *
psi_j_potential_plus
-
beta *
psi_i_potential_plus *
psi_j_field_plus
+
(scratch.penalty/h) *
psi_i_potential_plus *
psi_j_potential_plus
) *
scratch.fe_face_values.JxW(q);
}
}
}
}
else {
double h = std::min(interior_cell->diameter(),
exterior_cell->diameter());
scratch.fe_subface_values.reinit(interior_cell, interior_face_index, interior_subface_index);
scratch.fe_neighbour_face_values.reinit(exterior_cell, exterior_face_index);
const unsigned int n_face_points = scratch.fe_subface_values.n_quadrature_points;
const unsigned int dofs_this_cell = scratch.fe_subface_values.dofs_per_cell;
const unsigned int dofs_neighbor_cell = scratch.fe_neighbour_face_values.dofs_per_cell;
data.cell_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.cell_rhs.reinit(dofs_this_cell);
data.local_dof_indices.resize(dofs_this_cell);
data.local_neighbor_dof_indices.resize(dofs_neighbor_cell);
interior_cell->get_dof_indices(data.local_dof_indices);
exterior_cell->get_dof_indices(data.local_neighbor_dof_indices);
data.vi_ui_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.vi_ue_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.ve_ui_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.ve_ue_matrix.reinit(dofs_this_cell, dofs_this_cell);
data.cell_matrix = 0;
data.cell_rhs = 0;
data.vi_ui_matrix = 0;
data.vi_ue_matrix = 0;
data.ve_ui_matrix = 0;
data.ve_ue_matrix = 0;
const FEValuesExtractors::Vector VectorField(0);
const FEValuesExtractors::Scalar Potential(dim);
Point<dim> beta;
for (int i=0; i<dim; i++)
beta(i) = 1.0;
beta /= sqrt(beta.square() );
for (unsigned int q=0; q<n_face_points; q++)
{
for (unsigned int i=0; i<dofs_this_cell; i++)
{
const Tensor<1,dim> psi_i_field_minus =
scratch.fe_subface_values[VectorField].value(i,q);
const double psi_i_potential_minus =
scratch.fe_subface_values[Potential].value(i,q);
for (unsigned int j=0; j<dofs_this_cell; j++)
{
const Tensor<1,dim> psi_j_field_minus =
scratch.fe_subface_values[VectorField].value(j,q);
const double psi_j_potential_minus =
scratch.fe_subface_values[Potential].value(j,q);
data.vi_ui_matrix(i,j) += (0.5 * (
psi_i_field_minus *
scratch.fe_subface_values.normal_vector(q) *
psi_j_potential_minus
+
psi_i_potential_minus *
scratch.fe_subface_values.normal_vector(q) *
psi_j_field_minus )
+
beta *
psi_i_field_minus *
psi_j_potential_minus
-
beta *
psi_i_potential_minus *
psi_j_field_minus
+
(scratch.penalty/h) *
psi_i_potential_minus *
psi_j_potential_minus
) *
scratch.fe_subface_values.JxW(q);
}
for (unsigned int j=0; j<dofs_neighbor_cell; j++)
{
const Tensor<1,dim> psi_j_field_plus =
scratch.fe_neighbour_face_values[VectorField].value(j,q);
const double psi_j_potential_plus =
scratch.fe_neighbour_face_values[Potential].value(j,q);
data.vi_ue_matrix(i,j) += ( 0.5 * (
psi_i_field_minus *
scratch.fe_subface_values.normal_vector(q) *
psi_j_potential_plus
+
psi_i_potential_minus *
scratch.fe_subface_values.normal_vector(q) *
psi_j_field_plus )
-
beta *
psi_i_field_minus *
psi_j_potential_plus
+
beta *
psi_i_potential_minus *
psi_j_field_plus
-
(scratch.penalty/h) *
psi_i_potential_minus *
psi_j_potential_plus
) *
scratch.fe_subface_values.JxW(q);
}
}
for (unsigned int i=0; i<dofs_neighbor_cell; i++)
{
const Tensor<1,dim> psi_i_field_plus =
scratch.fe_neighbour_face_values[VectorField].value(i,q);
const double psi_i_potential_plus =
scratch.fe_neighbour_face_values[Potential].value(i,q);
for (unsigned int j=0; j<dofs_this_cell; j++)
{
const Tensor<1,dim> psi_j_field_minus =
scratch.fe_subface_values[VectorField].value(j,q);
const double psi_j_potential_minus =
scratch.fe_subface_values[Potential].value(j,q);
data.ve_ui_matrix(i,j) += (-0.5 * (
psi_i_field_plus *
scratch.fe_subface_values.normal_vector(q) *
psi_j_potential_minus
+
psi_i_potential_plus *
scratch.fe_subface_values.normal_vector(q) *
psi_j_field_minus)
-
beta *
psi_i_field_plus *
psi_j_potential_minus
+
beta *
psi_i_potential_plus *
psi_j_field_minus
-
(scratch.penalty/h) *
psi_i_potential_plus *
psi_j_potential_minus
) *
scratch.fe_subface_values.JxW(q);
}
for (unsigned int j=0; j<dofs_neighbor_cell; j++)
{
const Tensor<1,dim> psi_j_field_plus =
scratch.fe_neighbour_face_values[VectorField].value(j,q);
const double psi_j_potential_plus =
scratch.fe_neighbour_face_values[Potential].value(j,q);
data.ve_ue_matrix(i,j) += (-0.5 * (
psi_i_field_plus *
scratch.fe_face_values.normal_vector(q) *
psi_j_potential_plus
+
psi_i_potential_plus *
scratch.fe_face_values.normal_vector(q) *
psi_j_field_plus )
+
beta *
psi_i_field_plus *
psi_j_potential_plus
-
beta *
psi_i_potential_plus *
psi_j_field_plus
+
(scratch.penalty/h) *
psi_i_potential_plus *
psi_j_potential_plus
) *
scratch.fe_subface_values.JxW(q);
}
}
}
}
}
template <int dim>
void LDGPoissonProblem<dim>::local_assemble_boundaries(const typename DoFHandler<dim>::active_cell_iterator &cell,
const unsigned int face_index,
AssemblyScratchData &scratch,
AssemblyCopyData &data)
{
scratch.fe_face_values.reinit(cell, face_index);
if(cell->face(face_index)->boundary_id() == Dirichlet)
{
//Dirichlet
const unsigned int dofs_per_cell = scratch.fe_face_values.dofs_per_cell;
const unsigned int n_q_points = scratch.fe_face_values.n_quadrature_points;
const double h = cell->diameter();
data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.vi_ui_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.vi_ue_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.ve_ui_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.ve_ue_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.cell_rhs.reinit(dofs_per_cell);
data.local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(data.local_dof_indices);
data.cell_rhs = 0;
data.cell_matrix = 0;
data.vi_ui_matrix = 0;
data.vi_ue_matrix = 0;
data.ve_ui_matrix = 0;
data.ve_ue_matrix = 0;
const FEValuesExtractors::Vector VectorField(0);
const FEValuesExtractors::Scalar Potential(dim);
std::vector<double> Dirichlet_bc_values(n_q_points);
Dirichlet_bc_function.value_list(scratch.fe_face_values.get_quadrature_points(),
Dirichlet_bc_values);
for (unsigned int q=0; q < n_q_points; ++q)
{
for (unsigned int i=0; i < dofs_per_cell; ++i)
{
const Tensor<1, dim> psi_i_field = scratch.fe_face_values[VectorField].value(i,q);
const double psi_i_potential = scratch.fe_face_values[Potential].value(i,q);
for (unsigned int j=0; j<dofs_per_cell; j++)
{
const Tensor<1, dim> psi_j_field = scratch.fe_face_values[VectorField].value(j,q);
const double psi_j_potential = scratch.fe_face_values[Potential].value(j,q);
data.cell_matrix(i,j) += psi_i_potential * (
scratch.fe_face_values.normal_vector(q) *
psi_j_field
+
(scratch.penalty/h) *
psi_j_potential) *
scratch.fe_face_values.JxW(q);
}
data.cell_rhs(i) += (-1.0 * psi_i_field *
scratch.fe_face_values.normal_vector(q)
+
(scratch.penalty/h) *
psi_i_potential) *
Dirichlet_bc_values[q] *
scratch.fe_face_values.JxW(q);
}
}
}
else {
//Neumann
const unsigned int dofs_per_cell = scratch.fe_face_values.dofs_per_cell;
const unsigned int n_q_points = scratch.fe_face_values.n_quadrature_points;
data.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
data.cell_rhs.reinit(dofs_per_cell);
data.local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(data.local_dof_indices);
data.cell_rhs = 0;
data.cell_matrix = 0;
const FEValuesExtractors::Vector VectorField(0);
const FEValuesExtractors::Scalar Potential(dim);
std::vector<double > Neumann_bc_values(n_q_points);
for (unsigned int q = 0; q < n_q_points; ++q)
{
for (unsigned int i = 0; i<dofs_per_cell; ++i)
{
const Tensor<1, dim> psi_i_field = scratch.fe_face_values[VectorField].value(i,q);
const double psi_i_potential = scratch.fe_face_values[Potential].value(i,q);
for (unsigned int j = 0; j < dofs_per_cell; ++j)
{
const double psi_j_potential = scratch.fe_face_values[Potential].value(j,q);
data.cell_matrix(i,j) += psi_i_field *
scratch.fe_face_values.normal_vector(q) *
psi_j_potential *
scratch.fe_face_values.JxW(q);
}
data.cell_rhs(i) += -psi_i_potential *
Neumann_bc_values[q] *
scratch.fe_face_values.JxW(q);
}
}
}
}and for copying everything together
template <int dim>
void LDGPoissonProblem<dim>::copy_local_to_global_system(const AssemblyCopyData &data)
{
if(data.cell_matrix.m() == data.local_dof_indices.size())
constraints.distribute_local_to_global(data.cell_matrix,
data.cell_rhs,
data.local_dof_indices,
system_matrix,
system_rhs);
if(data.vi_ui_matrix.m() == data.local_dof_indices.size())
constraints.distribute_local_to_global(data.vi_ui_matrix,
data.local_dof_indices,
system_matrix);
if(data.vi_ue_matrix.m() == data.local_dof_indices.size() && data.vi_ue_matrix.n() == data.local_neighbor_dof_indices.size())
constraints.distribute_local_to_global(data.vi_ue_matrix,
data.local_dof_indices,
data.local_neighbor_dof_indices,
system_matrix);
if(data.ve_ui_matrix.m() == data.local_neighbor_dof_indices.size() && data.ve_ui_matrix.n() == data.local_dof_indices.size())
constraints.distribute_local_to_global(data.ve_ui_matrix,
data.local_neighbor_dof_indices,
data.local_dof_indices,
system_matrix);
if(data.ve_ue_matrix.m() == data.local_neighbor_dof_indices.size())
constraints.distribute_local_to_global(data.ve_ue_matrix,
data.local_neighbor_dof_indices,
system_matrix);
}
void LDGPoissonProblem<dim>::copy_local_to_global_system(const AssemblyCopyData &data)
{
if(data.cell_matrix.m() == data.local_dof_indices.size())
constraints.distribute_local_to_global(data.cell_matrix,
data.cell_rhs,
data.local_dof_indices,
system_matrix,
system_rhs);
if(data.vi_ui_matrix.m() == data.local_dof_indices.size())
constraints.distribute_local_to_global(data.vi_ui_matrix,
data.local_dof_indices,
system_matrix);
if(data.vi_ue_matrix.m() == data.local_dof_indices.size() && data.vi_ue_matrix.n() == data.local_neighbor_dof_indices.size())
constraints.distribute_local_to_global(data.vi_ue_matrix,
data.local_dof_indices,
data.local_neighbor_dof_indices,
system_matrix);
if(data.ve_ui_matrix.m() == data.local_neighbor_dof_indices.size() && data.ve_ui_matrix.n() == data.local_dof_indices.size())
constraints.distribute_local_to_global(data.ve_ui_matrix,
data.local_neighbor_dof_indices,
data.local_dof_indices,
system_matrix);
if(data.ve_ue_matrix.m() == data.local_neighbor_dof_indices.size())
constraints.distribute_local_to_global(data.ve_ue_matrix,
data.local_neighbor_dof_indices,
system_matrix);
}
If I omit the if()-guards, I will get the following error:
--------------------------------------------------------
An error occurred in line <3912> of file </home/roland/Downloads/git-files/dealii/include/deal.II/lac/affine_constraints.templates.h> in function
void dealii::AffineConstraints<number>::distribute_local_to_global(const dealii::FullMatrix<number>&, const std::vector<unsigned int>&, const dealii::AffineConstraints<number>&, const std::vector<unsigned int>&, MatrixType&) const [with MatrixType = dealii::TrilinosWrappers::SparseMatrix; number = double]
The violated condition was:
static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(local_matrix.n()), decltype(col_indices.size())>::type)>::type>(local_matrix.n()) == static_cast<typename ::dealii::internal::argument_type<void( typename std::common_type<decltype(local_matrix.n()), decltype(col_indices.size())>::type)>::type>(col_indices.size())
Additional information:
Dimension 12 not equal to 0.if I do not omit it, I get:
An error occurred in line <841> of file </home/roland/Downloads/git-files/dealii/source/lac/trilinos_solver.cc> in function
void dealii::TrilinosWrappers::SolverDirect::do_solve()
The violated condition was:
ierr == 0
Additional information:
An error with error number -22 occurred while calling a Trilinos functionThus I assume I either make a mistake when initializing the fluxmatrices, or when assembling the whole system. I tried debugging by printing the system matrix after assembling the system, but calling
std::ofstream out("system_matrix_meshworker.vtu");
matrix_out.build_patches(system_matrix, "system_matrix_meshworker");
matrix_out.write_vtu(out);led to std::bad_alloc, thus I am currently stuck in terms of how I could debug further. Is my general approach correct, or did I do something fundamentally wrong?
Thanks!
Maxi Miller
Mar 5, 2019, 5:53:35 AM3/5/19
to deal.II User Group
As initial debug approach, I used the function local_assemble_system instead of assemble_cell_terms, which worked (left the other functions alone). Then I tried to assemble just the boundary terms, by adding local_assemble_boundaries, uncommenting assemble_Dirichlet_boundary_terms() and assemble_Neumann_boundary_terms() and adding MeshWorker::AssembleFlags::assemble_boundary_faces to the flags in the mesh_loop. Here it failed, but I was not able to find out why yet.
Reply all
Reply to author
Forward
0 new messages