Assemble the system matrix and rhs vector.
212 const DofMap & dof_map = system.get_dof_map();
214 FEType fe_type = dof_map.variable_type(0);
216 std::unique_ptr<FEBase> fe (FEBase::build(
dim, fe_type));
217 std::unique_ptr<FEBase> fe_elem_face (FEBase::build(
dim, fe_type));
218 std::unique_ptr<FEBase> fe_neighbor_face (FEBase::build(
dim, fe_type));
223 fe->attach_quadrature_rule (&qrule);
224 fe_elem_face->attach_quadrature_rule (&qface);
225 fe_neighbor_face->attach_quadrature_rule (&qface);
227 const std::vector<Real> & JxW = fe->get_JxW();
228 const std::vector<std::vector<Real>> & phi = fe->get_phi();
229 const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
231 const std::vector<Real> & JxW_face = fe_elem_face->get_JxW();
233 const std::vector<Point> & qface_points = fe_elem_face->get_xyz();
235 const std::vector<std::vector<Real>> & phi_face = fe_elem_face->get_phi();
236 const std::vector<std::vector<Real>> & phi_neighbor_face = fe_neighbor_face->get_phi();
246 std::vector<dof_id_type> dof_indices;
249 for (
const auto & elem :
mesh.active_local_element_ptr_range())
251 dof_map.dof_indices (elem, dof_indices);
252 const unsigned int n_dofs = dof_indices.size();
256 Ke.
resize (n_dofs, n_dofs);
260 for (
unsigned int qp=0; qp<qrule.n_points(); qp++)
261 for (
unsigned int i=0; i<n_dofs; i++)
262 for (
unsigned int j=0; j<n_dofs; j++)
263 Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
267 for (
auto side : elem->side_index_range())
268 if (elem->neighbor_ptr(side) ==
nullptr)
272 fe_elem_face->reinit(elem, side);
274 for (
unsigned int qp=0; qp<qface.n_points(); qp++)
275 for (std::size_t i=0; i<phi.size(); i++)
276 Fe(i) += JxW_face[qp] * phi_face[i][qp];
284 for (
auto side : elem->side_index_range())
285 if (elem->neighbor_ptr(side) ==
nullptr)
290 fe_elem_face->reinit(elem, side);
292 ElementSideMap::const_iterator ltu_it =
293 lower_to_upper.find(std::make_pair(elem, side));
296 const Elem * neighbor = ltu_it->second;
298 std::vector<Point> qface_neighbor_points;
299 FEMap::inverse_map (elem->dim(), neighbor,
301 qface_neighbor_points);
302 fe_neighbor_face->reinit(neighbor, &qface_neighbor_points);
304 std::vector<dof_id_type> neighbor_dof_indices;
305 dof_map.dof_indices (neighbor, neighbor_dof_indices);
306 const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
308 Kne.
resize (n_neighbor_dofs, n_dofs);
309 Ken.
resize (n_dofs, n_neighbor_dofs);
310 Kee.
resize (n_dofs, n_dofs);
311 Knn.
resize (n_neighbor_dofs, n_neighbor_dofs);
314 for (
unsigned int qp=0; qp<qface.n_points(); qp++)
315 for (
unsigned int i=0; i<n_dofs; i++)
316 for (
unsigned int j=0; j<n_dofs; j++)
317 Kee(i,j) -= JxW_face[qp] * (1./R)*(phi_face[i][qp] * phi_face[j][qp]);
320 for (
unsigned int qp=0; qp<qface.n_points(); qp++)
321 for (
unsigned int i=0; i<n_dofs; i++)
322 for (
unsigned int j=0; j<n_neighbor_dofs; j++)
323 Ken(i,j) += JxW_face[qp] * (1./R)*(phi_face[i][qp] * phi_neighbor_face[j][qp]);
326 for (
unsigned int qp=0; qp<qface.n_points(); qp++)
327 for (
unsigned int i=0; i<n_neighbor_dofs; i++)
328 for (
unsigned int j=0; j<n_neighbor_dofs; j++)
329 Knn(i,j) -= JxW_face[qp] * (1./R)*(phi_neighbor_face[i][qp] * phi_neighbor_face[j][qp]);
332 for (
unsigned int qp=0; qp<qface.n_points(); qp++)
333 for (
unsigned int i=0; i<n_neighbor_dofs; i++)
334 for (
unsigned int j=0; j<n_dofs; j++)
335 Kne(i,j) += JxW_face[qp] * (1./R)*(phi_neighbor_face[i][qp] * phi_face[j][qp]);
337 matrix.
add_matrix(Kne, neighbor_dof_indices, dof_indices);
338 matrix.
add_matrix(Ken, dof_indices, neighbor_dof_indices);
345 dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
348 system.rhs->add_vector (Fe, dof_indices);
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
void resize(const unsigned int n)
Resize the vector.
This is the base class from which all geometric element types are derived.
Order default_quadrature_order() const
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
This class handles the numbering of degrees of freedom on a mesh.
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
const T & get(std::string_view) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const MeshBase & get_mesh() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Parameters parameters
Data structure holding arbitrary parameters.