libMesh
dg_fem_context.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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(libmesh_nullptr),
32  _neighbor_dof_indices_var(sys.n_vars()),
33  _dg_terms_active(false)
34 {
35  libmesh_experimental();
36 
37  unsigned int nv = sys.n_vars();
38  libmesh_assert (nv);
39 
40  _neighbor_subresiduals.reserve(nv);
41  _elem_elem_subjacobians.resize(nv);
42  _elem_neighbor_subjacobians.resize(nv);
43  _neighbor_elem_subjacobians.resize(nv);
45 
46  for (unsigned int i=0; i != nv; ++i)
47  {
49  _elem_elem_subjacobians[i].reserve(nv);
50  _elem_neighbor_subjacobians[i].reserve(nv);
51  _neighbor_elem_subjacobians[i].reserve(nv);
52  _neighbor_neighbor_subjacobians[i].reserve(nv);
53 
54  for (unsigned int j=0; j != nv; ++j)
55  {
56  _elem_elem_subjacobians[i].push_back
58  _elem_neighbor_subjacobians[i].push_back
60  _neighbor_elem_subjacobians[i].push_back
64  }
65  }
66 
67  _neighbor_side_fe_var.resize(nv);
68  for (unsigned int i=0; i != nv; ++i)
69  {
70  FEType fe_type = sys.variable_type(i);
71 
72  if (_neighbor_side_fe[fe_type] == libmesh_nullptr)
73  _neighbor_side_fe[fe_type] = FEAbstract::build(this->_dim, fe_type).release();
74 
76  }
77 }
78 
80 {
81 
82  for (std::size_t i=0; i != _neighbor_subresiduals.size(); ++i)
83  {
84  delete _neighbor_subresiduals[i];
85 
86  for (std::size_t j=0; j != _elem_elem_subjacobians[i].size(); ++j)
87  {
88  delete _elem_elem_subjacobians[i][j];
89  delete _elem_neighbor_subjacobians[i][j];
90  delete _neighbor_elem_subjacobians[i][j];
92  }
93  }
94 
95  // Delete the FE objects
96  for (std::map<FEType, FEAbstract *>::iterator i = _neighbor_side_fe.begin();
97  i != _neighbor_side_fe.end(); ++i)
98  delete i->second;
99  _neighbor_side_fe.clear();
100 }
101 
103 {
105 
106  // By default we assume that the DG terms are inactive
107  // They are only active if neighbor_side_fe_reinit is called
108  _dg_terms_active = false;
109 }
110 
112 {
113  // Call this *after* side_fe_reinit
114 
115  // Initialize all the neighbor side FE objects based on inverse mapping
116  // the quadrature points on the current side
117  std::vector<Point> qface_side_points;
118  std::vector<Point> qface_neighbor_points;
119  std::map<FEType, FEAbstract *>::iterator local_fe_end = _neighbor_side_fe.end();
120  for (std::map<FEType, FEAbstract *>::iterator i = _neighbor_side_fe.begin();
121  i != local_fe_end; ++i)
122  {
123  FEType neighbor_side_fe_type = i->first;
124  FEAbstract * side_fe = _side_fe[this->get_dim()][neighbor_side_fe_type];
125  qface_side_points = side_fe->get_xyz();
126 
128  neighbor_side_fe_type,
129  &get_neighbor(),
130  qface_side_points,
131  qface_neighbor_points);
132 
133  i->second->reinit(&get_neighbor(), &qface_neighbor_points);
134  }
135 
136  // Set boolean flag to indicate that the DG terms are active on this element
137  _dg_terms_active = true;
138 
139  // Also, initialize data required for DG assembly on the current side,
140  // analogously to FEMContext::pre_fe_reinit
141 
142  // Initialize the per-element data for elem.
144 
145  const unsigned int n_dofs = cast_int<unsigned int>
146  (this->get_dof_indices().size());
147  const unsigned int n_neighbor_dofs = cast_int<unsigned int>
148  (_neighbor_dof_indices.size());
149 
150  // These resize calls also zero out the residual and jacobian
151  _neighbor_residual.resize(n_neighbor_dofs);
152  _elem_elem_jacobian.resize(n_dofs, n_dofs);
153  _elem_neighbor_jacobian.resize(n_dofs, n_neighbor_dofs);
154  _neighbor_elem_jacobian.resize(n_neighbor_dofs, n_dofs);
155  _neighbor_neighbor_jacobian.resize(n_neighbor_dofs, n_neighbor_dofs);
156 
157  // Initialize the per-variable data for elem.
158  {
159  unsigned int sub_dofs = 0;
160  for (unsigned int i=0; i != get_system().n_vars(); ++i)
161  {
163 
164  const unsigned int n_dofs_var = cast_int<unsigned int>
165  (_neighbor_dof_indices_var[i].size());
166 
167  _neighbor_subresiduals[i]->reposition
168  (sub_dofs, n_dofs_var);
169 
170  for (unsigned int j=0; j != i; ++j)
171  {
172  const unsigned int n_dofs_var_j =
173  cast_int<unsigned int>
174  (this->get_dof_indices(j).size());
175 
176  _elem_elem_subjacobians[i][j]->reposition
177  (sub_dofs, _neighbor_subresiduals[j]->i_off(),
178  n_dofs_var, n_dofs_var_j);
179  _elem_elem_subjacobians[j][i]->reposition
180  (_neighbor_subresiduals[j]->i_off(), sub_dofs,
181  n_dofs_var_j, n_dofs_var);
182 
183  _elem_neighbor_subjacobians[i][j]->reposition
184  (sub_dofs, _neighbor_subresiduals[j]->i_off(),
185  n_dofs_var, n_dofs_var_j);
186  _elem_neighbor_subjacobians[j][i]->reposition
187  (_neighbor_subresiduals[j]->i_off(), sub_dofs,
188  n_dofs_var_j, n_dofs_var);
189 
190  _neighbor_elem_subjacobians[i][j]->reposition
191  (sub_dofs, _neighbor_subresiduals[j]->i_off(),
192  n_dofs_var, n_dofs_var_j);
193  _neighbor_elem_subjacobians[j][i]->reposition
194  (_neighbor_subresiduals[j]->i_off(), sub_dofs,
195  n_dofs_var_j, n_dofs_var);
196 
197  _neighbor_neighbor_subjacobians[i][j]->reposition
198  (sub_dofs, _neighbor_subresiduals[j]->i_off(),
199  n_dofs_var, n_dofs_var_j);
200  _neighbor_neighbor_subjacobians[j][i]->reposition
201  (_neighbor_subresiduals[j]->i_off(), sub_dofs,
202  n_dofs_var_j, n_dofs_var);
203  }
204  _elem_elem_subjacobians[i][i]->reposition
205  (sub_dofs, sub_dofs,
206  n_dofs_var,
207  n_dofs_var);
208  _elem_neighbor_subjacobians[i][i]->reposition
209  (sub_dofs, sub_dofs,
210  n_dofs_var,
211  n_dofs_var);
212  _neighbor_elem_subjacobians[i][i]->reposition
213  (sub_dofs, sub_dofs,
214  n_dofs_var,
215  n_dofs_var);
216  _neighbor_neighbor_subjacobians[i][i]->reposition
217  (sub_dofs, sub_dofs,
218  n_dofs_var,
219  n_dofs_var);
220  sub_dofs += n_dofs_var;
221  }
222  libmesh_assert_equal_to (sub_dofs, n_dofs);
223  }
224 
225 }
226 
227 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
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.
ImplicitSystem & sys
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
std::vector< std::vector< dof_id_type > > _neighbor_dof_indices_var
std::map< FEType, FEAbstract * > _neighbor_side_fe
Finite element objects for each variable&#39;s sides on the neighbor element.
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
const class libmesh_nullptr_t libmesh_nullptr
std::vector< FEAbstract * > _neighbor_side_fe_var
Pointers to the same finite element objects on the neighbor element, but indexed by variable number...
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
const FEType & variable_type(const unsigned int i) const
Definition: system.h:2164
const Elem & get_neighbor() const
Accessor for current neighbor Elem object for assembling DG terms.
libmesh_assert(j)
Defines a dense subvector for use in finite element computations.
unsigned char get_dim() const
Accessor for cached mesh dimension.
Definition: fem_context.h:899
const DofMap & get_dof_map() const
Definition: system.h:2030
std::vector< std::vector< DenseSubMatrix< Number > * > > _elem_elem_subjacobians
virtual void side_fe_reinit()
Reinitializes side FE objects on the current geometric element.
Definition: fem_context.C:1405
std::vector< std::vector< DenseSubMatrix< Number > * > > _neighbor_neighbor_subjacobians
std::vector< std::vector< DenseSubMatrix< Number > * > > _elem_neighbor_subjacobians
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
Defines a dense submatrix for use in Finite Element-type computations.
DenseMatrix< Number > _neighbor_elem_jacobian
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:569
DenseMatrix< Number > _elem_neighbor_jacobian
const System & get_system() const
Accessor for associated system.
Definition: diff_context.h:104
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:1122
std::vector< std::map< FEType, FEAbstract * > > _side_fe
Definition: fem_context.h:1094
std::vector< std::vector< DenseSubMatrix< Number > * > > _neighbor_elem_subjacobians
DGFEMContext(const System &sys)
Constructor.
DenseMatrix< Number > _neighbor_neighbor_jacobian
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
std::vector< DenseSubVector< Number > * > _neighbor_subresiduals
Element residual subvectors and Jacobian submatrices.
unsigned int n_vars() const
Definition: system.h:2086
virtual ~DGFEMContext()
Destructor.
static UniquePtr< FEAbstract > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
Definition: fe_abstract.C:44
This class forms the foundation from which generic finite elements may be derived.
Definition: fe_abstract.h:93
virtual void side_fe_reinit() libmesh_override
Override side_fe_reinit to set a boolean flag so that by default DG terms are assumed to be inactive...
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
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:1917