Franco Milicchio
the issue with MeshWorker is that when it was implemented, we didn't write enough documentation/tutorial programs/examples [...]. Development then also stopped, and so [...] nobody today knows any more whether they can or can not be done with MeshWorker. The only person who still knows how this all works unfortunately
also doesn't regularly participate in the help forum any more :-(
[...] you can of course ask questions about how this or that can be done, but you're likely not going to get a lot of helpful answers because nobody knows any more how this works. It's really an unfortunate situation we got ourselves into with this small corner of the library.
// Integrate local cell matrix
void integrate_cell_term(MeshWorker::DoFInfo<dim> &dinfo, MeshWorker::IntegrationInfo<dim> &info)
{
loopc++;
// FE Values
const FEValuesBase<dim> &fe_v = info.fe_values();
// Local Lame values, here are stupid constants
std::vector<double> lambda_values(fe_v.n_quadrature_points), mu_values(fe_v.n_quadrature_points);
// Local element matrix
FullMatrix<double> &local_matrix = dinfo.matrix(0).matrix;
// Local element RHS vector
BlockVector<double> &local_vector = dinfo.vector(0);
// Jacobian
const std::vector<double> &JxW = fe_v.get_JxW_values ();
// Local Lame values
lambda.value_list(fe_v.get_quadrature_points(), lambda_values);
mu.value_list(fe_v.get_quadrature_points(), mu_values);
// For each DOF of a cell (i-th)
for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
{
const unsigned int component_i = info.fe_values().get_fe().system_to_component_index(i).first;
// For each DOF of a cell (j-th)
for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
{
const unsigned int component_j = info.fe_values().get_fe().system_to_component_index(j).first;
// Compute the (i, j)-th component for LHS via quadrature
for (unsigned int q_point = 0; q_point < fe_v.n_quadrature_points; ++q_point)
{
local_matrix(i,j) += ((fe_v.shape_grad(i, q_point)[component_i] *
fe_v.shape_grad(j, q_point)[component_j] *
lambda_values[q_point])
+
(fe_v.shape_grad(i, q_point)[component_j] *
fe_v.shape_grad(j, q_point)[component_i] *
mu_values[q_point])
+
(
(i == j) ?
(fe_v.shape_grad(component_i, q_point) *
fe_v.shape_grad(component_j, q_point) *
mu_values[q_point])
: 0)
) * JxW[q_point];
} // Integral for bilinear form
} // j
// Computed values for the RHS
std::vector<Vector<double>> rhs_values(fe_v.n_quadrature_points, Vector<double>(dim));
// Computed values for the RHS
std::vector<Vector<double>> rhs_face_values(fe_v.n_quadrature_points, Vector<double>(dim));
// RHS values
right_hand_side.vector_value_list(fe_v.get_quadrature_points(), rhs_values);
// Compute the i-th component for RHS via quadrature (body forces)
for(unsigned int q_point = 0; q_point < fe_v.n_quadrature_points; ++q_point)
{
local_vector(i) += fe_v.shape_value(i, q_point) *
rhs_values[q_point](component_i) *
JxW[q_point];
} // Integral for linear form
} // i
// Computed values for the RHS on the faces
std::vector<Vector<double>> rhs_face_values(fe_v.n_quadrature_points, Vector<double>(dim));
// Neumann (boundary forces on top boundary == id 2)
for (unsigned int face_number = 0; face_number < GeometryInfo<dim>::faces_per_cell; ++face_number)
{
loopf++;
if (dinfo.cell->face(face_number)->at_boundary() && (dinfo.cell->face(face_number)->boundary_id() == 2))
{
// FE Values on face (???)
const FEValuesBase<dim> &fe_face_values = info.fe_values();
// Get face FE values for this face and cell
fe_face_values.reinit(dinfo.cell, face_number);
// Neumann force vector on the face
right_face_side.vector_value_list(fe_face_values.get_quadrature_points(), rhs_face_values);
for(unsigned int f_q_point = 0; f_q_point < fe_v.n_quadrature_points; ++f_q_point)
{
// Add contributions to all DOFs
for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
{
const unsigned int component_i = fe_v.system_to_component_index(i).first;
local_vector(i) += (fe_face_values.shape_value(i, f_q_point) *
rhs_face_values[f_q_point][component_i]) * fe_face_values.JxW(f_q_point);
}
} // Integral for linear form
}
} // Boundary forces
}
Daniel Arndt
First of all, I have read this question, and I read that MeshWorker is sort of abandoned:the issue with MeshWorker is that when it was implemented, we didn't write enough documentation/tutorial programs/examples [...]. Development then also stopped, and so [...] nobody today knows any more whether they can or can not be done with MeshWorker. The only person who still knows how this all works unfortunately
also doesn't regularly participate in the help forum any more :-(
[...] you can of course ask questions about how this or that can be done, but you're likely not going to get a lot of helpful answers because nobody knows any more how this works. It's really an unfortunate situation we got ourselves into with this small corner of the library.So my first question is, should I avoid using this class and implement parallel loops by hand (via TBB or other means)?
The problem really is the documentation, but with working code you should be able to understand the concepts.
Then, my doubts about the use of MeshWorkers. I cannot really grasp it even after reading tutorials 12, 16, and 39.What is the purpose of the cell, boundary, and face functions? The cell is called exactly on each n-dimensional cell, I've checked that number, what about the other two? Boundary is called, as I read, for each n-dimensional cell on the boundary (numbers again match). Face function is called on each internal (n-1)-dimensional face (not on the boundary), as far as I understand. Is this correct?
In a typical DG formulation your bilinear form consists of a cell term, a boundary term and face term. These relate to the cell, boundary and face functions you are talking about. You should be able to see this in step-39.
When assigning a Neumann condition on faces, for instance in elasticity apply a load, where is the correct place? I have tried to get the code into the cell function, but I didn't find a way to get the faces, boundary markers, and the face values. My code is at the bottom, if you need it.
Neumann boundary conditions are normally considered in a weak sense and incorporated as boundary term into the bilinear form. So this should be in the boundary method of your MatrixIntegrator.
There is an elasticity example in amandus. Maybe this is useful.
Is there an analogous of the FullMatrix for the right hand side vector? The only thing similar to a vector in the DoFInfo is a BlockVector...
If you initialize the "MeshWorker::DoFInfo" object with a "BlockInfo" object, you also have access to individual blocks. The structure then resembles the structure of your FiniteElement.
Franco Milicchio
On Monday, March 6, 2017 at 7:49:53 PM UTC+1, Daniel Arndt wrote:
So my first question is, should I avoid using this class and implement parallel loops by hand (via TBB or other means)?
"amandus"[1] is in fact based on MeshWorker. If you are trying to do something similar to what is in the examples there, you have a good chance to be successful.
The problem really is the documentation, but with working code you should be able to understand the concepts.
Then, my doubts about the use of MeshWorkers. I cannot really grasp it even after reading tutorials 12, 16, and 39.What is the purpose of the cell, boundary, and face functions? The cell is called exactly on each n-dimensional cell, I've checked that number, what about the other two? Boundary is called, as I read, for each n-dimensional cell on the boundary (numbers again match). Face function is called on each internal (n-1)-dimensional face (not on the boundary), as far as I understand. Is this correct?
In MeshWorker, you are basically only writing the assembly loops for your bilinear forms. You don't have to care about how to handle hanging nodes and the like.
In a typical DG formulation your bilinear form consists of a cell term, a boundary term and face term. These relate to the cell, boundary and face functions you are talking about. You should be able to see this in step-39.
When assigning a Neumann condition on faces, for instance in elasticity apply a load, where is the correct place? I have tried to get the code into the cell function, but I didn't find a way to get the faces, boundary markers, and the face values. My code is at the bottom, if you need it.
In general, your code fragment looks good, but the right-hand side should be handled in a different integrator object similar to what is done in step-39.
Neumann boundary conditions are normally considered in a weak sense and incorporated as boundary term into the bilinear form. So this should be in the boundary method of your MatrixIntegrator.
There is an elasticity example in amandus. Maybe this is useful.
Is there an analogous of the FullMatrix for the right hand side vector? The only thing similar to a vector in the DoFInfo is a BlockVector...
What exactly do you mean? In general, all vectors are considered to be block vectors with only one block that you can access via "dinfo.vector(0).block(0);".
If you initialize the "MeshWorker::DoFInfo" object with a "BlockInfo" object, you also have access to individual blocks. The structure then resembles the structure of your FiniteElement.
Franco Milicchio
static void rhs_integrate_boundary_term (MeshWorker::DoFInfo<dim> &dinfo, MeshWorker::IntegrationInfo<dim> &info)
{
if (dumpall)
{
std::lock_guard<std::mutex> guard(lock);
std::cout << "> RHS boundary " << dinfo.cell->id() << " face " << dinfo.face_number << std::endl;
}
loopb++;
// FE Values
const FEValuesBase<dim> &fe_v = info.fe_values();
// Local RHS vector
Vector<double> &local_vector = dinfo.vector(0).block(0);
// Local values
std::vector<Vector<double>> boundary_values(fe_v.n_quadrature_points, Vector<double>(dim));
// Compute values for this boundary cell
right_face_side.vector_value_list(fe_v.get_quadrature_points(), boundary_values);
for (unsigned k = 0; k < fe_v.n_quadrature_points; ++k)
for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
{
const unsigned int component_i = info.fe_values().get_fe().system_to_component_index(i).first;
local_vector(i) += (fe_v.shape_value(i, k) * boundary_values[k][component_i]) * fe_v.JxW(k);
}
}
Franco Milicchio
Daniel Arndt
on the first glance your code looks good. Can you reduce the size of your problem to a single cell?
Is there still a difference then?
Best,
Daniel
Franco Milicchio
Franco Milicchio
Jean-Paul Pelteret
((i == j) ?
(fe_v.shape_grad(component_i, q_point) *
fe_v.shape_grad(component_j, q_point) *
mu_values[q_point]) :
0)
((component_i == component_j) ?
(fe_v.shape_grad(i, q_point) *
fe_v.shape_grad(j, q_point) *
mu_values[q_point]) :
0)