v0.14.0
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
22 auto fe_ptr = getPtrFE();
23 auto fe_type = getFEType();
24 // Set data structure to store base
25 if (dataL2->dataOnEntities[MBVERTEX].size() != 1) {
26 dataL2->dataOnEntities[MBVERTEX].clear();
27 dataL2->dataOnEntities[MBVERTEX].push_back(
29 }
30 if (dataL2->dataOnEntities[fe_type].size() != 1) {
31 dataL2->dataOnEntities[fe_type].clear();
32 dataL2->dataOnEntities[fe_type].push_back(new EntitiesFieldData::EntData());
33 }
34
35 auto &vertex_data = dataL2->dataOnEntities[MBVERTEX][0];
36 vertex_data.getNSharedPtr(NOBASE) =
37 fe_ptr->getEntData(H1, MBVERTEX, 0).getNSharedPtr(NOBASE);
38 vertex_data.getDiffNSharedPtr(NOBASE) =
39 fe_ptr->getEntData(H1, MBVERTEX, 0).getDiffNSharedPtr(NOBASE);
40
41 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
42 ent_data.getSense() = 1;
43 ent_data.getSpace() = L2;
44 ent_data.getBase() = base;
45 ent_data.getOrder() = get_order();
46
47 CHKERR fe_ptr->getElementPolynomialBase()->getValue(
48 getGaussPts(), boost::make_shared<EntPolynomialBaseCtx>(
49 *dataL2, static_cast<FieldSpace>(L2),
50 static_cast<FieldApproximationBase>(base), NOBASE));
52}
53
56 auto fe_type = getFEType();
57 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
58 auto &base_funcions = ent_data.getN(base);
59 const auto nb = base_funcions.size2();
60
61 if (nb) {
62
63 const auto nb_integration_pts = getGaussPts().size2();
64
65#ifndef NDEBUG
66 if (!baseMassPtr)
67 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
68 "Mass matrix is null pointer");
69#endif
70
71 auto &nN = *baseMassPtr;
72 nN.resize(nb, nb, false);
73 nN.clear();
74
76 // get base function gradient on rows
77 auto t_row_base = ent_data.getFTensor0N();
78 // loop over integration points
79 for (int gg = 0; gg != nb_integration_pts; ++gg) {
80 // take into account Jacobian
81 const double alpha = t_w;
82
83 for (int rr = 0; rr != nb; ++rr) {
84
85 // loop over rows base functions
86 auto a_mat_ptr = &nN(rr, 0);
87 // get column base functions gradient at gauss point gg
88 auto t_col_base = ent_data.getFTensor0N(gg, 0);
89 // loop over columns
90 for (int cc = 0; cc <= rr; ++cc) {
91 // calculate element of local matrix
92 *a_mat_ptr += alpha * (t_row_base * t_col_base);
93 ++t_col_base;
94 ++a_mat_ptr;
95 }
96
97 ++t_row_base;
98 }
99 ++t_w; // move to another integration weight
100 }
101
103 }
105}
106
111
112 if (sPace != L2) {
113 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
114 "Space should be set to L2");
115 }
116
117 CHKERR calculateBase(
118 [this]() { return std::max(0, getPtrFE()->getMaxDataOrder() - 1); });
119 CHKERR calculateMass();
120
122}
123
125 boost::shared_ptr<EntitiesFieldData> data_l2,
126 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
127 : OpSetInvJacSpaceForFaceImpl<2, 1>(NOSPACE, inv_jac_ptr), dataL2(data_l2) {
128}
129
134
135 auto apply_transform = [&](MatrixDouble &diff_n) {
136 return applyTransform<2, 2, 2, 2>(diff_n);
137 };
138
139 const auto fe_type = getFEType();
140 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
141 CHKERR apply_transform(ent_data.getDiffN());
142
144}
145
147 int derivative, boost::shared_ptr<MatrixDouble> base_mass_ptr,
148 boost::shared_ptr<EntitiesFieldData> data_l2,
149 const FieldApproximationBase b, const FieldSpace s, int verb, Sev sev)
150 : OpBaseDerivativesBase(base_mass_ptr, data_l2, b, s, verb, sev),
151 calcBaseDerivative(derivative) {
152 if (s != H1)
153 doEntities[MBVERTEX] = false;
154}
155
156template <int BASE_DIM, int SPACE_DIM>
159 EntitiesFieldData::EntData &ent_data) {
161
162 const int nb_gauss_pts = data.getN(base).size1();
163 const int nb_approx_bases = data.getN(base).size2() / BASE_DIM;
164 const int nb_derivatives =
165 BASE_DIM * std::pow(SPACE_DIM, calcBaseDerivative - 1);
166
167 const int nb_prj_bases = ent_data.getN().size2();
168
169 auto &n_diff_shared_ptr = data.getNSharedPtr(
170 base, static_cast<BaseDerivatives>(calcBaseDerivative));
171 if (!n_diff_shared_ptr)
172 n_diff_shared_ptr = boost::make_shared<MatrixDouble>();
173
174 auto &nex_diff_base = *(n_diff_shared_ptr);
175 const int next_nb_derivatives = BASE_DIM * pow(SPACE_DIM, calcBaseDerivative);
176 nex_diff_base.resize(nb_gauss_pts, nb_approx_bases * next_nb_derivatives,
177 false);
178 nex_diff_base.clear();
179
181 auto next_diffs_ptr = &*nex_diff_base.data().begin();
182 auto t_next_diff = getFTensor1FromPtr<SPACE_DIM>(next_diffs_ptr);
183
184 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
185
186 auto ptr = &*nF.data().begin();
187
188 for (auto r = 0; r != nb_approx_bases * nb_derivatives; ++r) {
189
190 auto l2_diff_base = ent_data.getFTensor1DiffN<SPACE_DIM>(base, gg, 0);
191 for (int rr = 0; rr != nb_prj_bases; ++rr) {
192 t_next_diff(i) += l2_diff_base(i) * (*ptr);
193 ++l2_diff_base;
194 ++ptr;
195 }
196
197 ++t_next_diff;
198 }
199 }
200
202}
203
204template <int BASE_DIM>
209
210 auto &approx_base = data.getN(base);
211 const auto nb_approx_bases = approx_base.size2() / BASE_DIM;
212
213 if (nb_approx_bases) {
214
215 const auto fe_type = getFEType();
216 const auto nb_integration_pts = approx_base.size1();
217
218 const auto space_dim =
219 data.getDiffN(base).size2() / (BASE_DIM * nb_approx_bases);
220 auto &diff_approx_base = *(data.getNSharedPtr(
221 base, static_cast<BaseDerivatives>(calcBaseDerivative - 1)));
222 int nb_derivatives = BASE_DIM * pow(space_dim, calcBaseDerivative - 1);
223
224 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
225 const int nb_prj_bases = ent_data.getN().size2();
226
227#ifndef NDEBUG
228 if (!baseMassPtr)
229 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
230 "Mass matrix is null pointer");
231#endif
232 auto &nN = *baseMassPtr;
233
234#ifndef NDEBUG
235 if (diff_approx_base.size2() != nb_approx_bases * nb_derivatives) {
236 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
237 "Number of derivatives and basses do not match");
238 }
239 if (ent_data.getN().size1() != nb_integration_pts) {
240 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
241 "Number of integration points is not consistent");
242 }
243 if (nN.size2() != nb_prj_bases) {
244 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
245 "Number of base functions and size of mass matrix does not math");
246 }
247#endif
248
249 nF.resize(nb_approx_bases * nb_derivatives, nb_prj_bases, false);
250 nF.clear();
251
252 auto t_w = getFTensor0IntegrationWeight();
253 // get base function gradient on rows
254
255 auto diff_base_ptr = &*diff_approx_base.data().begin();
256 // loop over integration points
257 for (int gg = 0; gg != nb_integration_pts; ++gg) {
258 // take into account Jacobian
259 const double alpha = t_w;
260
261 for (int r = 0; r != nb_approx_bases * nb_derivatives; ++r) {
262
263 // Rows are base functions
264 auto t_base = ent_data.getFTensor0N(base, gg, 0);
265 for (int rr = 0; rr != nb_prj_bases; ++rr) {
266 nF(r, rr) += alpha * (t_base * (*diff_base_ptr));
267 ++t_base;
268 }
269
270 ++diff_base_ptr;
271 }
272
273 ++t_w; // move to another integration weight
274 }
275
276 for (auto r = 0; r != nb_approx_bases * nb_derivatives; ++r) {
277 ublas::matrix_row<MatrixDouble> mc(nF, r);
278 cholesky_solve(nN, mc, ublas::lower());
279 }
280
281 if (space_dim == 3)
282 CHKERR setBaseImpl<BASE_DIM, 3>(data, ent_data);
283 else if (space_dim == 2)
284 CHKERR setBaseImpl<BASE_DIM, 2>(data, ent_data);
285 // else if (space_dim == 1)
286 // CHKERR setBaseImpl<BASE_DIM, 1>(data, ent_data);
287 else
288 SETERRQ1(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE,
289 "Space dim can be only 1,2,3 but is %d", space_dim);
290 }
291
293}
294
298 return doWorkImpl<1>(side, type, data);
299}
300
304 return doWorkImpl<3>(side, type, data);
305}
306
307} // 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
constexpr int BASE_DIM
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
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)
EntityType getFEType() const
Get dimension of finite element.
auto getFTensor0IntegrationWeight()
Get integration weights.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
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)
boost::shared_ptr< EntitiesFieldData > dataL2
MoFEMErrorCode calculateBase(GetOrderFun get_order)
boost::shared_ptr< MatrixDouble > baseMassPtr
Transform local reference derivatives of shape functions to global derivatives.