v0.13.2
Loading...
Searching...
No Matches
BaseDerivativesDataOperators.cpp
Go to the documentation of this file.
1/** \file BaseDerivativesDataOperators.cpp
2
3
4*/
5
6
7
8#include <cholesky.hpp>
9
10namespace MoFEM {
11
13 boost::shared_ptr<MatrixDouble> base_mass_ptr,
14 boost::shared_ptr<EntitiesFieldData> data_l2,
15 const FieldApproximationBase b, const FieldSpace s, int verb, Sev sev)
16 : ForcesAndSourcesCore::UserDataOperator(s, OPSPACE), base(b),
17 verbosity(verb), severityLevel(sev), baseMassPtr(base_mass_ptr),
18 dataL2(data_l2) {}
19
24
25 if (sPace != L2) {
26 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
27 "Space should be set to L2");
28 }
29
30 auto fe_type = getFEType();
31
32 auto calculate_base = [&]() {
34 auto fe_ptr = getPtrFE();
35 // Set data structure to store base
36 dataL2->dataOnEntities[fe_type].clear();
37 dataL2->dataOnEntities[MBVERTEX].push_back(
39 dataL2->dataOnEntities[fe_type].push_back(new EntitiesFieldData::EntData());
40
41 auto &vertex_data = dataL2->dataOnEntities[MBVERTEX][0];
42 vertex_data.getNSharedPtr(NOBASE) =
43 fe_ptr->getEntData(H1, MBVERTEX, 0).getNSharedPtr(NOBASE);
44 vertex_data.getDiffNSharedPtr(NOBASE) =
45 fe_ptr->getEntData(H1, MBVERTEX, 0).getDiffNSharedPtr(NOBASE);
46
47 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
48 ent_data.getSense() = 1;
49 ent_data.getBase() = base;
50 ent_data.getOrder() = std::max(0, fe_ptr->getMaxDataOrder() - 1);
51
52 CHKERR fe_ptr->getElementPolynomialBase()->getValue(
53 getGaussPts(), boost::make_shared<EntPolynomialBaseCtx>(
54 *dataL2, static_cast<FieldSpace>(L2),
55 static_cast<FieldApproximationBase>(base), NOBASE));
57 };
58
59 CHKERR calculate_base();
60
61 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
62 auto &base_funcions = ent_data.getN(base);
63 const auto nb = base_funcions.size2();
64
65 if (nb) {
66
67 const auto nb_integration_pts = getGaussPts().size2();
68
69#ifndef NDEBUG
70 if (!baseMassPtr)
71 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
72 "Mass matrix is null pointer");
73#endif
74
75 auto &nN = *baseMassPtr;
76 nN.resize(nb, nb, false);
77 nN.clear();
78
79 auto t_w = getFTensor0IntegrationWeight();
80 // get base function gradient on rows
81 auto t_row_base = ent_data.getFTensor0N();
82 // loop over integration points
83 for (int gg = 0; gg != nb_integration_pts; ++gg) {
84 // take into account Jacobian
85 const double alpha = t_w;
86
87 for (int rr = 0; rr != nb; ++rr) {
88
89 // loop over rows base functions
90 auto a_mat_ptr = &nN(rr, 0);
91 // get column base functions gradient at gauss point gg
92 auto t_col_base = ent_data.getFTensor0N(gg, 0);
93 // loop over columns
94 for (int cc = 0; cc <= rr; ++cc) {
95 // calculate element of local matrix
96 *a_mat_ptr += alpha * (t_row_base * t_col_base);
97 ++t_col_base;
98 ++a_mat_ptr;
99 }
100
101 ++t_row_base;
102 }
103 ++t_w; // move to another integration weight
104 }
105
107 }
108
110}
111
113 boost::shared_ptr<EntitiesFieldData> data_l2,
114 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
115 : OpSetInvJacSpaceForFaceImpl<2, 1>(NOSPACE, inv_jac_ptr), dataL2(data_l2) {
116}
117
122
123 auto apply_transform = [&](MatrixDouble &diff_n) {
124 return applyTransform<2, 2, 2, 2>(diff_n);
125 };
126
127 const auto fe_type = getFEType();
128 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
129 CHKERR apply_transform(ent_data.getDiffN());
130
132}
133
135 int derivative, boost::shared_ptr<MatrixDouble> base_mass_ptr,
136 boost::shared_ptr<EntitiesFieldData> data_l2,
137 const FieldApproximationBase b, const FieldSpace s, int verb, Sev sev)
138 : OpBaseDerivativesBase(base_mass_ptr, data_l2, b, s, verb, sev),
139 calcBaseDerivative(derivative) {
140 if (s != H1)
141 doEntities[MBVERTEX] = false;
142}
143
144template <int BASE_DIM, int SPACE_DIM>
147 EntitiesFieldData::EntData &ent_data) {
149
150 const int nb_gauss_pts = data.getN(base).size1();
151 const int nb_approx_bases = data.getN(base).size2() / BASE_DIM;
152 const int nb_derivatives =
153 BASE_DIM * std::pow(SPACE_DIM, calcBaseDerivative - 1);
154
155 const int nb_prj_bases = ent_data.getN().size2();
156
157 auto &n_diff_shared_ptr = data.getNSharedPtr(
158 base, static_cast<BaseDerivatives>(calcBaseDerivative));
159 if (!n_diff_shared_ptr)
160 n_diff_shared_ptr = boost::make_shared<MatrixDouble>();
161
162 auto &nex_diff_base = *(n_diff_shared_ptr);
163 const int next_nb_derivatives = BASE_DIM * pow(SPACE_DIM, calcBaseDerivative);
164 nex_diff_base.resize(nb_gauss_pts, nb_approx_bases * next_nb_derivatives,
165 false);
166 nex_diff_base.clear();
167
169 auto next_diffs_ptr = &*nex_diff_base.data().begin();
170 auto t_next_diff = getFTensor1FromPtr<SPACE_DIM>(next_diffs_ptr);
171
172 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
173
174 auto ptr = &*nF.data().begin();
175
176 for (auto r = 0; r != nb_approx_bases * nb_derivatives; ++r) {
177
178 auto l2_diff_base = ent_data.getFTensor1DiffN<SPACE_DIM>(base, gg, 0);
179 for (int rr = 0; rr != nb_prj_bases; ++rr) {
180 t_next_diff(i) += l2_diff_base(i) * (*ptr);
181 ++l2_diff_base;
182 ++ptr;
183 }
184
185 ++t_next_diff;
186 }
187 }
188
190}
191
192template <int BASE_DIM>
197
198 auto &approx_base = data.getN(base);
199 const auto nb_approx_bases = approx_base.size2() / BASE_DIM;
200
201 if (nb_approx_bases) {
202
203 const auto fe_type = getFEType();
204 const auto nb_integration_pts = approx_base.size1();
205
206 const auto space_dim =
207 data.getDiffN(base).size2() / (BASE_DIM * nb_approx_bases);
208 auto &diff_approx_base = *(data.getNSharedPtr(
209 base, static_cast<BaseDerivatives>(calcBaseDerivative - 1)));
210 int nb_derivatives = BASE_DIM * pow(space_dim, calcBaseDerivative - 1);
211
212 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
213 const int nb_prj_bases = ent_data.getN().size2();
214
215#ifndef NDEBUG
216 if (!baseMassPtr)
217 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
218 "Mass matrix is null pointer");
219#endif
220 auto &nN = *baseMassPtr;
221
222#ifndef NDEBUG
223 if (diff_approx_base.size2() != nb_approx_bases * nb_derivatives) {
224 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
225 "Number of derivatives and basses do not match");
226 }
227 if (ent_data.getN().size1() != nb_integration_pts) {
228 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
229 "Number of integration points is not consistent");
230 }
231 if (nN.size2() != nb_prj_bases) {
232 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
233 "Number of base functions and size of mass matrix does not math");
234 }
235#endif
236
237 nF.resize(nb_approx_bases * nb_derivatives, nb_prj_bases, false);
238 nF.clear();
239
240 auto t_w = getFTensor0IntegrationWeight();
241 // get base function gradient on rows
242
243 auto diff_base_ptr = &*diff_approx_base.data().begin();
244 // loop over integration points
245 for (int gg = 0; gg != nb_integration_pts; ++gg) {
246 // take into account Jacobian
247 const double alpha = t_w;
248
249 for (int r = 0; r != nb_approx_bases * nb_derivatives; ++r) {
250
251 // Rows are base functions
252 auto t_base = ent_data.getFTensor0N(base, gg, 0);
253 for (int rr = 0; rr != nb_prj_bases; ++rr) {
254 nF(r, rr) += alpha * (t_base * (*diff_base_ptr));
255 ++t_base;
256 }
257
258 ++diff_base_ptr;
259 }
260
261 ++t_w; // move to another integration weight
262 }
263
264 for (auto r = 0; r != nb_approx_bases * nb_derivatives; ++r) {
265 ublas::matrix_row<MatrixDouble> mc(nF, r);
266 cholesky_solve(nN, mc, ublas::lower());
267 }
268
269 if (space_dim == 3)
270 CHKERR setBaseImpl<BASE_DIM, 3>(data, ent_data);
271 else if (space_dim == 2)
272 CHKERR setBaseImpl<BASE_DIM, 2>(data, ent_data);
273 // else if (space_dim == 1)
274 // CHKERR setBaseImpl<BASE_DIM, 1>(data, ent_data);
275 else
276 SETERRQ1(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
277 "Space dim can be only 1,2,3 but is %d", space_dim);
278 }
279
281}
282
286 return doWorkImpl<1>(side, type, data);
287}
288
292 return doWorkImpl<3>(side, type, data);
293}
294
295} // namespace MoFEM
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr int SPACE_DIM
cholesky decomposition
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ NOBASE
Definition: definitions.h:59
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ 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_IMPOSSIBLE_CASE
Definition: definitions.h:35
@ 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
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
constexpr int BASE_DIM
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
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....
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives direvatie)
structure to get information form mofem into EntitiesFieldData
OpBaseDerivativesBase(boost::shared_ptr< MatrixDouble > base_mass_ptr, boost::shared_ptr< EntitiesFieldData > data_l2, const FieldApproximationBase b, const FieldSpace s, int verb=QUIET, Sev sev=Sev::verbose)
Transform local reference derivatives of shape functions to global derivatives.