v0.14.0
Loading...
Searching...
No Matches
PoissonDiscontinousGalerkin.hpp
Go to the documentation of this file.
1/**
2 * \file PoissonDiscontinousGalerkin.hpp
3 * \example PoissonDiscontinousGalerkin.hpp
4 *
5 * Example of implementation for discontinuous Galerkin.
6 */
7
8
9
10// Define name if it has not been defined yet
11#ifndef __POISSON2DISGALERKIN_HPP__
12#define __POISSON2DISGALERKIN_HPP__
13
14// Namespace that contains necessary UDOs, will be included in the main program
16
17// Declare FTensor index for 2D problem
19
21// data for skeleton computation
22std::array<VectorInt, 2>
23 indicesRowSideMap; ///< indices on rows for left hand-side
24std::array<VectorInt, 2>
25 indicesColSideMap; ///< indices on columns for left hand-side
26std::array<MatrixDouble, 2> rowBaseSideMap; // base functions on rows
27std::array<MatrixDouble, 2> colBaseSideMap; // base function on columns
28std::array<MatrixDouble, 2> rowDiffBaseSideMap; // derivative of base functions
29std::array<MatrixDouble, 2> colDiffBaseSideMap; // derivative of base functions
30std::array<double, 2> areaMap; // area/volume of elements on the side
31std::array<int, 2> senseMap; // orientaton of local element edge/face in respect
32 // to global orientation of edge/face
33
34/**
35 * @brief Operator tp collect data from elements on the side of Edge/Face
36 *
37 */
39
40 OpCalculateSideData(std::string field_name, std::string col_field_name)
41 : FaceSideOp(field_name, col_field_name, FaceSideOp::OPROWCOL) {
42
43 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
44
45 for (auto t = moab::CN::TypeDimensionMap[SPACE_DIM].first;
46 t <= moab::CN::TypeDimensionMap[SPACE_DIM].second; ++t)
47 doEntities[t] = true;
48 }
49
50 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
51 EntityType col_type, EntData &row_data,
52 EntData &col_data) {
54
55 // Note: THat for L2 base data rows, and columns are the same, so operator
56 // above can be simpler operator for the right hand side, and data can be
57 // stored only for rows, since for columns data are the same. However for
58 // complex multi-physics problems that not necessary would be a case. For
59 // generality, we keep generic case.
60
61 if ((CN::Dimension(row_type) == SPACE_DIM) &&
62 (CN::Dimension(col_type) == SPACE_DIM)) {
63 const auto nb_in_loop = getFEMethod()->nInTheLoop;
64 indicesRowSideMap[nb_in_loop] = row_data.getIndices();
65 indicesColSideMap[nb_in_loop] = col_data.getIndices();
66 rowBaseSideMap[nb_in_loop] = row_data.getN();
67 colBaseSideMap[nb_in_loop] = col_data.getN();
68 rowDiffBaseSideMap[nb_in_loop] = row_data.getDiffN();
69 colDiffBaseSideMap[nb_in_loop] = col_data.getDiffN();
70 areaMap[nb_in_loop] = getMeasure();
71 senseMap[nb_in_loop] = getSkeletonSense();
72 if (!nb_in_loop) {
73 indicesRowSideMap[1].clear();
74 indicesColSideMap[1].clear();
75 rowBaseSideMap[1].clear();
76 colBaseSideMap[1].clear();
77 rowDiffBaseSideMap[1].clear();
78 colDiffBaseSideMap[1].clear();
79 areaMap[1] = 0;
80 senseMap[1] = 0;
81 }
82 } else {
83 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Should not happen");
84 }
85
87 }
88};
89
90template <typename T> inline auto get_ntensor(T &base_mat) {
92 &*base_mat.data().begin());
93};
94
95template <typename T> inline auto get_ntensor(T &base_mat, int gg, int bb) {
96 double *ptr = &base_mat(gg, bb);
98};
99
100template <typename T> inline auto get_diff_ntensor(T &base_mat) {
101 double *ptr = &*base_mat.data().begin();
102 return getFTensor1FromPtr<2>(ptr);
103};
104
105template <typename T>
106inline auto get_diff_ntensor(T &base_mat, int gg, int bb) {
107 double *ptr = &base_mat(gg, 2 * bb);
108 return getFTensor1FromPtr<2>(ptr);
109};
110
111/**
112 * @brief Operator the left hand side matrix
113 *
114 */
116public:
117
118 /**
119 * @brief Construct a new OpL2LhsPenalty
120 *
121 * @param side_fe_ptr pointer to FE to evaluate side elements
122 */
123 OpL2LhsPenalty(boost::shared_ptr<FaceSideEle> side_fe_ptr)
125
126 MoFEMErrorCode doWork(int side, EntityType type,
127 EntitiesFieldData::EntData &data) {
129
130 // Collect data from side domian elements
131 CHKERR loopSideFaces("dFE", *sideFEPtr);
132 const auto in_the_loop =
133 sideFEPtr->nInTheLoop; // return number of elements on the side
134
135#ifndef NDEBUG
136 const std::array<std::string, 2> ele_type_name = {"BOUNDARY", "SKELETON"};
137 MOFEM_LOG("SELF", Sev::noisy)
138 << "OpL2LhsPenalty inTheLoop " << ele_type_name[in_the_loop];
139 MOFEM_LOG("SELF", Sev::noisy)
140 << "OpL2LhsPenalty sense " << senseMap[0] << " " << senseMap[1];
141#endif
142
143 // calculate penalty
144 const double s = getMeasure() / (areaMap[0] + areaMap[1]);
145 const double p = penalty * s;
146
147 // get normal of the face or edge
148 auto t_normal = getFTensor1Normal();
149 t_normal.normalize();
150
151 // get number of integration points on face
152 const size_t nb_integration_pts = getGaussPts().size2();
153
154 // beta paramter is zero, when penalty method is used, also, takes value 1,
155 // when boundary edge/face is evaluated, and 2 when is skeleton edge/face.
156 const double beta = static_cast<double>(nitsche) / (in_the_loop + 1);
157
158 // iterate the sides rows
159 for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
160
161 // gent number of DOFs on the right side.
162 const auto nb_rows = indicesRowSideMap[s0].size();
163
164 if (nb_rows) {
165
166 // get orientation of the local element edge
167 const auto sense_row = senseMap[s0];
168
169 // iterate the side cols
170 const auto nb_row_base_functions = rowBaseSideMap[s0].size2();
171 for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
172
173 // get orientation of the edge
174 const auto sense_col = senseMap[s1];
175
176 // number of dofs, for homogenous approximation this value is
177 // constant.
178 const auto nb_cols = indicesColSideMap[s1].size();
179
180 // resize local element matrix
181 locMat.resize(nb_rows, nb_cols, false);
182 locMat.clear();
183
184 // get base functions, and integration weights
185 auto t_row_base = get_ntensor(rowBaseSideMap[s0]);
186 auto t_diff_row_base = get_diff_ntensor(rowDiffBaseSideMap[s0]);
187 auto t_w = getFTensor0IntegrationWeight();
188
189 // iterate integration points on face/edge
190 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
191
192 // t_w is integration weight, and measure is element area, or
193 // volume, depending if problem is in 2d/3d.
194 const double alpha = getMeasure() * t_w;
195 auto t_mat = locMat.data().begin();
196
197 // iterate rows
198 size_t rr = 0;
199 for (; rr != nb_rows; ++rr) {
200
201 // calculate tetting function
203 t_vn_plus(i) = beta * (phi * t_diff_row_base(i) / p);
205 t_vn(i) = t_row_base * t_normal(i) * sense_row - t_vn_plus(i);
206
207 // get base functions on columns
208 auto t_col_base = get_ntensor(colBaseSideMap[s1], gg, 0);
209 auto t_diff_col_base =
211
212 // iterate columns
213 for (size_t cc = 0; cc != nb_cols; ++cc) {
214
215 // calculate variance of tested function
217 t_un(i) = -p * (t_col_base * t_normal(i) * sense_col -
218 beta * t_diff_col_base(i) / p);
219
220 // assemble matrix
221 *t_mat -= alpha * (t_vn(i) * t_un(i));
222 *t_mat -= alpha * (t_vn_plus(i) * (beta * t_diff_col_base(i)));
223
224 // move to next column base and element of matrix
225 ++t_col_base;
226 ++t_diff_col_base;
227 ++t_mat;
228 }
229
230 // move to next row base
231 ++t_row_base;
232 ++t_diff_row_base;
233 }
234
235 // this is obsolete for this particular example, we keep it for
236 // generality. in case of multi-physics problems different fields can
237 // chare hierarchical base but use different approx. order, so is
238 // possible to have more base functions than DOFs on element.
239 for (; rr < nb_row_base_functions; ++rr) {
240 ++t_row_base;
241 ++t_diff_row_base;
242 }
243
244 ++t_w;
245 }
246
247 // assemble system
248 CHKERR ::MatSetValues(getKSPB(), indicesRowSideMap[s0].size(),
249 &*indicesRowSideMap[s0].begin(),
250 indicesColSideMap[s1].size(),
251 &*indicesColSideMap[s1].begin(),
252 &*locMat.data().begin(), ADD_VALUES);
253
254 // stop of boundary element
255 if (!in_the_loop)
257 }
258 }
259 }
260
262 }
263
264private:
265 boost::shared_ptr<FaceSideEle>
266 sideFEPtr; ///< pointer to element to get data on edge/face sides
267 MatrixDouble locMat; ///< local operator matrix
268};
269
270/**
271 * @brief Operator to evaluate Dirichlet boundary conditions using DG
272 *
273 */
275public:
276 OpL2BoundaryRhs(boost::shared_ptr<FaceSideEle> side_fe_ptr, ScalarFun fun)
278 sourceFun(fun) {}
279
280 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
282
283 // get normal of the face or edge
284 CHKERR loopSideFaces("dFE", *sideFEPtr);
285 const auto in_the_loop =
286 sideFEPtr->nInTheLoop; // return number of elements on the side
287
288 // calculate penalty
289 const double s = getMeasure() / (areaMap[0]);
290 const double p = penalty * s;
291
292 // get normal of the face or edge
293 auto t_normal = getFTensor1Normal();
294 t_normal.normalize();
295
296 auto t_w = getFTensor0IntegrationWeight();
297
298 const size_t nb_rows = indicesRowSideMap[0].size();
299
300 if (!nb_rows)
302
303 // resize local operator vector
304 rhsF.resize(nb_rows, false);
305 rhsF.clear();
306
307 const size_t nb_integration_pts = getGaussPts().size2();
308 const size_t nb_row_base_functions = rowBaseSideMap[0].size2();
309
310 // shape funcs
311 auto t_row_base = get_ntensor(rowBaseSideMap[0]);
312 auto t_diff_row_base = get_diff_ntensor(rowDiffBaseSideMap[0]);
313 auto t_coords = getFTensor1CoordsAtGaussPts();
314
315 const auto sense_row = senseMap[0];
316 const double beta = static_cast<double>(nitsche) / (in_the_loop + 1);
317
318 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
319 const double alpha = getMeasure() * t_w;
320
321 const double source_val =
322 -p * sourceFun(t_coords(0), t_coords(1), t_coords(2));
323
324 auto t_f = rhsF.data().begin();
325
326 size_t rr = 0;
327 for (; rr != nb_rows; ++rr) {
328
330 t_vn_plus(i) = beta * (phi * t_diff_row_base(i) / p);
332 t_vn(i) = t_row_base * t_normal(i) * sense_row - t_vn_plus(i);
333
334 // assemble operator local vactor
335 *t_f -= alpha * t_vn(i) * (source_val * t_normal(i));
336
337 ++t_row_base;
338 ++t_diff_row_base;
339 ++t_f;
340 }
341
342 for (; rr < nb_row_base_functions; ++rr) {
343 ++t_row_base;
344 ++t_diff_row_base;
345 }
346
347 ++t_w;
348 ++t_coords;
349 }
350
351 // assemble local operator vector to global vector
352 CHKERR ::VecSetValues(getKSPf(), indicesRowSideMap[0].size(),
353 &*indicesRowSideMap[0].begin(), &*rhsF.data().begin(),
354 ADD_VALUES);
355
357 }
358
359private:
360 boost::shared_ptr<FaceSideEle>
361 sideFEPtr; ///< pointer to element to get data on edge/face sides
362 ScalarFun sourceFun; ///< pointer to function to evaluate value of function on boundary
363 VectorDouble rhsF; ///< vector to store local operator right hand side
364};
365
366}; // namespace Poisson2DiscontGalerkinOperators
367
368#endif //__POISSON2DISGALERKIN_HPP__
static Index< 'p', 3 > p
constexpr int SPACE_DIM
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ NOSPACE
Definition: definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
auto fun
Function to approximate.
static double penalty
static double phi
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
const double T
std::array< MatrixDouble, 2 > rowBaseSideMap
std::array< VectorInt, 2 > indicesColSideMap
indices on columns for left hand-side
FTensor::Index< 'i', SPACE_DIM > i
std::array< MatrixDouble, 2 > colDiffBaseSideMap
std::array< MatrixDouble, 2 > colBaseSideMap
std::array< VectorInt, 2 > indicesRowSideMap
indices on rows for left hand-side
std::array< MatrixDouble, 2 > rowDiffBaseSideMap
constexpr double t
plate stiffness
Definition: plate.cpp:59
static double nitsche
constexpr auto field_name
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Operator tp collect data from elements on the side of Edge/Face.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpCalculateSideData(std::string field_name, std::string col_field_name)
Operator to evaluate Dirichlet boundary conditions using DG.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
VectorDouble rhsF
vector to store local operator right hand side
OpL2BoundaryRhs(boost::shared_ptr< FaceSideEle > side_fe_ptr, ScalarFun fun)
ScalarFun sourceFun
pointer to function to evaluate value of function on boundary
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
OpL2LhsPenalty(boost::shared_ptr< FaceSideEle > side_fe_ptr)
Construct a new OpL2LhsPenalty.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides