Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 :
29 3025 : DGFEMContext::DGFEMContext (const System & sys)
30 : : FEMContext(sys),
31 2857 : _neighbor(nullptr),
32 2941 : _neighbor_dof_indices_var(sys.n_vars()),
33 3109 : _dg_terms_active(false)
34 : {
35 84 : unsigned int nv = sys.n_vars();
36 84 : libmesh_assert (nv);
37 :
38 3025 : _neighbor_subresiduals.reserve(nv);
39 3025 : _elem_elem_subjacobians.resize(nv);
40 3025 : _elem_neighbor_subjacobians.resize(nv);
41 3025 : _neighbor_elem_subjacobians.resize(nv);
42 3025 : _neighbor_neighbor_subjacobians.resize(nv);
43 :
44 8890 : for (unsigned int i=0; i != nv; ++i)
45 : {
46 5865 : _neighbor_subresiduals.emplace_back(std::make_unique<DenseSubVector<Number>>(_neighbor_residual));
47 6029 : _elem_elem_subjacobians[i].reserve(nv);
48 6029 : _elem_neighbor_subjacobians[i].reserve(nv);
49 6029 : _neighbor_elem_subjacobians[i].reserve(nv);
50 6029 : _neighbor_neighbor_subjacobians[i].reserve(nv);
51 :
52 20250 : for (unsigned int j=0; j != nv; ++j)
53 : {
54 14385 : _elem_elem_subjacobians[i].emplace_back(std::make_unique<DenseSubMatrix<Number>>(_elem_elem_jacobian));
55 14385 : _elem_neighbor_subjacobians[i].emplace_back(std::make_unique<DenseSubMatrix<Number>>(_elem_neighbor_jacobian));
56 14385 : _neighbor_elem_subjacobians[i].emplace_back(std::make_unique<DenseSubMatrix<Number>>(_neighbor_elem_jacobian));
57 14789 : _neighbor_neighbor_subjacobians[i].emplace_back(std::make_unique<DenseSubMatrix<Number>>(_neighbor_neighbor_jacobian));
58 : }
59 : }
60 :
61 3025 : _neighbor_side_fe_var.resize(nv);
62 8890 : for (unsigned int i=0; i != nv; ++i)
63 : {
64 5865 : FEType fe_type = sys.variable_type(i);
65 :
66 5865 : if (_neighbor_side_fe[fe_type] == nullptr)
67 5966 : _neighbor_side_fe[fe_type] = FEAbstract::build(this->_dim, fe_type);
68 :
69 5865 : _neighbor_side_fe_var[i] = _neighbor_side_fe[fe_type].get();
70 : }
71 3025 : }
72 :
73 14151 : DGFEMContext::~DGFEMContext() = default;
74 :
75 136060 : void DGFEMContext::side_fe_reinit()
76 : {
77 136060 : FEMContext::side_fe_reinit();
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 136060 : _dg_terms_active = false;
82 136060 : }
83 :
84 14400 : void DGFEMContext::neighbor_side_fe_reinit ()
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 2400 : std::vector<Point> qface_side_points;
91 2400 : std::vector<Point> qface_neighbor_points;
92 28800 : for (auto & [neighbor_side_fe_type, fe] : _neighbor_side_fe)
93 : {
94 15600 : FEAbstract * side_fe = _side_fe[this->get_dim()][neighbor_side_fe_type].get();
95 14400 : qface_side_points = side_fe->get_xyz();
96 :
97 14400 : FEMap::inverse_map (this->get_dim(), &get_neighbor(),
98 : qface_side_points, qface_neighbor_points);
99 :
100 14400 : 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 14400 : _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.
110 15600 : get_system().get_dof_map().dof_indices (&get_neighbor(), _neighbor_dof_indices);
111 :
112 : const unsigned int n_dofs = cast_int<unsigned int>
113 2400 : (this->get_dof_indices().size());
114 : const unsigned int n_neighbor_dofs = cast_int<unsigned int>
115 2400 : (_neighbor_dof_indices.size());
116 :
117 : // These resize calls also zero out the residual and jacobian
118 13200 : _neighbor_residual.resize(n_neighbor_dofs);
119 13200 : _elem_elem_jacobian.resize(n_dofs, n_dofs);
120 13200 : _elem_neighbor_jacobian.resize(n_dofs, n_neighbor_dofs);
121 13200 : _neighbor_elem_jacobian.resize(n_neighbor_dofs, n_dofs);
122 13200 : _neighbor_neighbor_jacobian.resize(n_neighbor_dofs, n_neighbor_dofs);
123 :
124 : // Initialize the per-variable data for elem.
125 : {
126 1200 : unsigned int sub_dofs = 0;
127 28800 : for (auto i : make_range(get_system().n_vars()))
128 : {
129 16800 : get_system().get_dof_map().dof_indices (&get_neighbor(), _neighbor_dof_indices_var[i], i);
130 :
131 : const unsigned int n_dofs_var = cast_int<unsigned int>
132 3600 : (_neighbor_dof_indices_var[i].size());
133 :
134 2400 : _neighbor_subresiduals[i]->reposition
135 1200 : (sub_dofs, n_dofs_var);
136 :
137 14400 : for (unsigned int j=0; j != i; ++j)
138 : {
139 : const unsigned int n_dofs_var_j =
140 : cast_int<unsigned int>
141 0 : (this->get_dof_indices(j).size());
142 :
143 0 : _elem_elem_subjacobians[i][j]->reposition
144 0 : (sub_dofs, _neighbor_subresiduals[j]->i_off(),
145 : n_dofs_var, n_dofs_var_j);
146 0 : _elem_elem_subjacobians[j][i]->reposition
147 0 : (_neighbor_subresiduals[j]->i_off(), sub_dofs,
148 : n_dofs_var_j, n_dofs_var);
149 :
150 0 : _elem_neighbor_subjacobians[i][j]->reposition
151 0 : (sub_dofs, _neighbor_subresiduals[j]->i_off(),
152 : n_dofs_var, n_dofs_var_j);
153 0 : _elem_neighbor_subjacobians[j][i]->reposition
154 0 : (_neighbor_subresiduals[j]->i_off(), sub_dofs,
155 : n_dofs_var_j, n_dofs_var);
156 :
157 0 : _neighbor_elem_subjacobians[i][j]->reposition
158 0 : (sub_dofs, _neighbor_subresiduals[j]->i_off(),
159 : n_dofs_var, n_dofs_var_j);
160 0 : _neighbor_elem_subjacobians[j][i]->reposition
161 0 : (_neighbor_subresiduals[j]->i_off(), sub_dofs,
162 : n_dofs_var_j, n_dofs_var);
163 :
164 0 : _neighbor_neighbor_subjacobians[i][j]->reposition
165 0 : (sub_dofs, _neighbor_subresiduals[j]->i_off(),
166 : n_dofs_var, n_dofs_var_j);
167 0 : _neighbor_neighbor_subjacobians[j][i]->reposition
168 0 : (_neighbor_subresiduals[j]->i_off(), sub_dofs,
169 : n_dofs_var_j, n_dofs_var);
170 : }
171 3600 : _elem_elem_subjacobians[i][i]->reposition
172 1200 : (sub_dofs, sub_dofs,
173 : n_dofs_var,
174 : n_dofs_var);
175 3600 : _elem_neighbor_subjacobians[i][i]->reposition
176 1200 : (sub_dofs, sub_dofs,
177 : n_dofs_var,
178 : n_dofs_var);
179 3600 : _neighbor_elem_subjacobians[i][i]->reposition
180 1200 : (sub_dofs, sub_dofs,
181 : n_dofs_var,
182 : n_dofs_var);
183 3600 : _neighbor_neighbor_subjacobians[i][i]->reposition
184 1200 : (sub_dofs, sub_dofs,
185 : n_dofs_var,
186 : n_dofs_var);
187 14400 : sub_dofs += n_dofs_var;
188 : }
189 1200 : libmesh_assert_equal_to (sub_dofs, n_dofs);
190 : }
191 :
192 14400 : }
193 :
194 : }
|
