libMesh
dg_fem_context.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "libmesh/dg_fem_context.h"
19 #include "libmesh/dof_map.h"
20 #include "libmesh/elem.h"
21 #include "libmesh/fe_base.h"
22 #include "libmesh/fe_interface.h"
23 #include "libmesh/quadrature.h"
24 #include "libmesh/system.h"
25 
26 namespace libMesh
27 {
28 
30  : FEMContext(sys),
31  _neighbor(nullptr),
32  _neighbor_dof_indices_var(sys.n_vars()),
33  _dg_terms_active(false)
34 {
35  unsigned int nv = sys.n_vars();
36  libmesh_assert (nv);
37 
38  _neighbor_subresiduals.reserve(nv);
39  _elem_elem_subjacobians.resize(nv);
40  _elem_neighbor_subjacobians.resize(nv);
41  _neighbor_elem_subjacobians.resize(nv);
43 
44  for (unsigned int i=0; i != nv; ++i)
45  {
46  _neighbor_subresiduals.emplace_back(std::make_unique<DenseSubVector<Number>>(_neighbor_residual));
47  _elem_elem_subjacobians[i].reserve(nv);
48  _elem_neighbor_subjacobians[i].reserve(nv);
49  _neighbor_elem_subjacobians[i].reserve(nv);
50  _neighbor_neighbor_subjacobians[i].reserve(nv);
51 
52  for (unsigned int j=0; j != nv; ++j)
53  {
54  _elem_elem_subjacobians[i].emplace_back(std::make_unique<DenseSubMatrix<Number>>(_elem_elem_jacobian));
58  }
59  }
60 
61  _neighbor_side_fe_var.resize(nv);
62  for (unsigned int i=0; i != nv; ++i)
63  {
64  FEType fe_type = sys.variable_type(i);
65 
66  if (_neighbor_side_fe[fe_type] == nullptr)
67  _neighbor_side_fe[fe_type] = FEAbstract::build(this->_dim, fe_type);
68 
69  _neighbor_side_fe_var[i] = _neighbor_side_fe[fe_type].get();
70  }
71 }
72 
73 DGFEMContext::~DGFEMContext() = default;
74 
76 {
78 
79  // By default we assume that the DG terms are inactive
80  // They are only active if neighbor_side_fe_reinit is called
81  _dg_terms_active = false;
82 }
83 
85 {
86  // Call this *after* side_fe_reinit
87 
88  // Initialize all the neighbor side FE objects based on inverse mapping
89  // the quadrature points on the current side
90  std::vector<Point> qface_side_points;
91  std::vector<Point> qface_neighbor_points;
92  for (auto & [neighbor_side_fe_type, fe] : _neighbor_side_fe)
93  {
94  FEAbstract * side_fe = _side_fe[this->get_dim()][neighbor_side_fe_type].get();
95  qface_side_points = side_fe->get_xyz();
96 
98  qface_side_points, qface_neighbor_points);
99 
100  fe->reinit(&get_neighbor(), &qface_neighbor_points);
101  }
102 
103  // Set boolean flag to indicate that the DG terms are active on this element
104  _dg_terms_active = true;
105 
106  // Also, initialize data required for DG assembly on the current side,
107  // analogously to FEMContext::pre_fe_reinit
108 
109  // Initialize the per-element data for elem.
111 
112  const unsigned int n_dofs = cast_int<unsigned int>
113  (this->get_dof_indices().size());
114  const unsigned int n_neighbor_dofs = cast_int<unsigned int>
115  (_neighbor_dof_indices.size());
116 
117  // These resize calls also zero out the residual and jacobian
118  _neighbor_residual.resize(n_neighbor_dofs);
119  _elem_elem_jacobian.resize(n_dofs, n_dofs);
120  _elem_neighbor_jacobian.resize(n_dofs, n_neighbor_dofs);
121  _neighbor_elem_jacobian.resize(n_neighbor_dofs, n_dofs);
122  _neighbor_neighbor_jacobian.resize(n_neighbor_dofs, n_neighbor_dofs);
123 
124  // Initialize the per-variable data for elem.
125  {
126  unsigned int sub_dofs = 0;
127  for (auto i : make_range(get_system().n_vars()))
128  {
130 
131  const unsigned int n_dofs_var = cast_int<unsigned int>
132  (_neighbor_dof_indices_var[i].size());
133 
134  _neighbor_subresiduals[i]->reposition
135  (sub_dofs, n_dofs_var);
136 
137  for (unsigned int j=0; j != i; ++j)
138  {
139  const unsigned int n_dofs_var_j =
140  cast_int<unsigned int>
141  (this->get_dof_indices(j).size());
142 
143  _elem_elem_subjacobians[i][j]->reposition
144  (sub_dofs, _neighbor_subresiduals[j]->i_off(),
145  n_dofs_var, n_dofs_var_j);
146  _elem_elem_subjacobians[j][i]->reposition
147  (_neighbor_subresiduals[j]->i_off(), sub_dofs,
148  n_dofs_var_j, n_dofs_var);
149 
150  _elem_neighbor_subjacobians[i][j]->reposition
151  (sub_dofs, _neighbor_subresiduals[j]->i_off(),
152  n_dofs_var, n_dofs_var_j);
153  _elem_neighbor_subjacobians[j][i]->reposition
154  (_neighbor_subresiduals[j]->i_off(), sub_dofs,
155  n_dofs_var_j, n_dofs_var);
156 
157  _neighbor_elem_subjacobians[i][j]->reposition
158  (sub_dofs, _neighbor_subresiduals[j]->i_off(),
159  n_dofs_var, n_dofs_var_j);
160  _neighbor_elem_subjacobians[j][i]->reposition
161  (_neighbor_subresiduals[j]->i_off(), sub_dofs,
162  n_dofs_var_j, n_dofs_var);
163 
164  _neighbor_neighbor_subjacobians[i][j]->reposition
165  (sub_dofs, _neighbor_subresiduals[j]->i_off(),
166  n_dofs_var, n_dofs_var_j);
167  _neighbor_neighbor_subjacobians[j][i]->reposition
168  (_neighbor_subresiduals[j]->i_off(), sub_dofs,
169  n_dofs_var_j, n_dofs_var);
170  }
171  _elem_elem_subjacobians[i][i]->reposition
172  (sub_dofs, sub_dofs,
173  n_dofs_var,
174  n_dofs_var);
175  _elem_neighbor_subjacobians[i][i]->reposition
176  (sub_dofs, sub_dofs,
177  n_dofs_var,
178  n_dofs_var);
179  _neighbor_elem_subjacobians[i][i]->reposition
180  (sub_dofs, sub_dofs,
181  n_dofs_var,
182  n_dofs_var);
183  _neighbor_neighbor_subjacobians[i][i]->reposition
184  (sub_dofs, sub_dofs,
185  n_dofs_var,
186  n_dofs_var);
187  sub_dofs += n_dofs_var;
188  }
189  libmesh_assert_equal_to (sub_dofs, n_dofs);
190  }
191 
192 }
193 
194 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
std::vector< dof_id_type > _neighbor_dof_indices
Global Degree of freedom index lists for the neighbor element.
DenseVector< Number > _neighbor_residual
Residual vector of the neighbor component.
DenseMatrix< Number > _elem_elem_jacobian
The DG Jacobian terms.
void neighbor_side_fe_reinit()
Initialize neighbor side data needed to assemble DG terms.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1992
std::vector< std::vector< std::unique_ptr< DenseSubMatrix< Number > > > > _elem_elem_subjacobians
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition: fe_map.C:1626
virtual void side_fe_reinit() override
Override side_fe_reinit to set a boolean flag so that by default DG terms are assumed to be inactive...
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:374
std::vector< std::vector< dof_id_type > > _neighbor_dof_indices_var
std::vector< FEAbstract * > _neighbor_side_fe_var
Pointers to the same finite element objects on the neighbor element, but indexed by variable number...
The libMesh namespace provides an interface to certain functionality in the library.
unsigned char get_dim() const
Accessor for cached mesh dimension.
Definition: fem_context.h:936
Defines a dense subvector for use in finite element computations.
virtual void side_fe_reinit()
Reinitializes side FE objects on the current geometric element.
Definition: fem_context.C:1490
const System & get_system() const
Accessor for associated system.
Definition: diff_context.h:106
unsigned int n_vars
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
Defines a dense submatrix for use in Finite Element-type computations.
DenseMatrix< Number > _neighbor_elem_jacobian
libmesh_assert(ctx)
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:62
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:363
DenseMatrix< Number > _elem_neighbor_jacobian
virtual_for_inffe const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:272
std::vector< std::vector< std::unique_ptr< DenseSubMatrix< Number > > > > _elem_neighbor_subjacobians
bool _dg_terms_active
Boolean flag to indicate whether or not the DG terms have been assembled and should be used in the gl...
unsigned char _dim
Cached dimension of largest dimension element in this mesh.
Definition: fem_context.h:1197
std::vector< std::map< FEType, std::unique_ptr< FEAbstract > > > _side_fe
Definition: fem_context.h:1169
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2427
std::vector< std::vector< std::unique_ptr< DenseSubMatrix< Number > > > > _neighbor_elem_subjacobians
DGFEMContext(const System &sys)
Constructor.
static std::unique_ptr< FEAbstract > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
Definition: fe_abstract.C:77
DenseMatrix< Number > _neighbor_neighbor_jacobian
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:134
std::vector< std::unique_ptr< DenseSubVector< Number > > > _neighbor_subresiduals
Element residual subvectors and Jacobian submatrices.
virtual ~DGFEMContext()
Destructor.
unsigned int n_vars() const
Definition: system.h:2349
This class forms the foundation from which generic finite elements may be derived.
Definition: fe_abstract.h:99
unsigned int n_vars() const
Number of variables in solution.
Definition: diff_context.h:100
std::map< FEType, std::unique_ptr< FEAbstract > > _neighbor_side_fe
Finite element objects for each variable&#39;s sides on the neighbor element.
const DofMap & get_dof_map() const
Definition: system.h:2293
const Elem & get_neighbor() const
Accessor for current neighbor Elem object for assembling DG terms.
std::vector< std::vector< std::unique_ptr< DenseSubMatrix< Number > > > > _neighbor_neighbor_subjacobians