dg_fem_context.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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(NULL),
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] == NULL )
73  {
74  _neighbor_side_fe[fe_type] = FEAbstract::build(dim, fe_type).release();
75  }
77  }
78 }
79 
81 {
82 
83  for (std::size_t i=0; i != _neighbor_subresiduals.size(); ++i)
84  {
85  delete _neighbor_subresiduals[i];
86 
87  for (std::size_t j=0; j != _elem_elem_subjacobians[i].size(); ++j)
88  {
89  delete _elem_elem_subjacobians[i][j];
90  delete _elem_neighbor_subjacobians[i][j];
91  delete _neighbor_elem_subjacobians[i][j];
93  }
94  }
95 
96  // Delete the FE objects
97  for (std::map<FEType, FEAbstract *>::iterator i = _neighbor_side_fe.begin();
98  i != _neighbor_side_fe.end(); ++i)
99  delete i->second;
100  _neighbor_side_fe.clear();
101 }
102 
104 {
106 
107  // By default we assume that the DG terms are inactive
108  // They are only active if neighbor_side_fe_reinit is called
109  _dg_terms_active = false;
110 }
111 
113 {
114  // Call this *after* side_fe_reinit
115 
116  // Initialize all the neighbor side FE objects based on inverse mapping
117  // the quadrature points on the current side
118  std::vector<Point> qface_side_points;
119  std::vector<Point> qface_neighbor_points;
120  std::map<FEType, FEAbstract *>::iterator local_fe_end = _neighbor_side_fe.end();
121  for (std::map<FEType, FEAbstract *>::iterator i = _neighbor_side_fe.begin();
122  i != local_fe_end; ++i)
123  {
124  FEType neighbor_side_fe_type = i->first;
125  FEAbstract* side_fe = _side_fe[neighbor_side_fe_type];
126  qface_side_points = side_fe->get_xyz();
127 
129  neighbor_side_fe_type,
130  &get_neighbor(),
131  qface_side_points,
132  qface_neighbor_points);
133 
134  i->second->reinit(&get_neighbor(), &qface_neighbor_points);
135  }
136 
137  // Set boolean flag to indicate that the DG terms are active on this element
138  _dg_terms_active = true;
139 
140  // Also, initialize data required for DG assembly on the current side,
141  // analogously to FEMContext::pre_fe_reinit
142 
143  // Initialize the per-element data for elem.
145 
146  const unsigned int n_dofs = dof_indices.size();
147  const unsigned int n_neighbor_dofs = libmesh_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 = libmesh_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  libmesh_cast_int<unsigned int>
174  (dof_indices_var[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 }
228 

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo