v0.16.3
Loading...
Searching...
No Matches
BiLinearFormsIntegratorsImpl.hpp
Go to the documentation of this file.
1/** \file BiLinearFormsIntegratorsImpl.hpp
2 * \brief Bilinear forms integrators (implementation)
3 * \ingroup mofem_form
4
5 \todo Some operators could be optimised. To do that, we need to write tests
6 and use Valgrind to profile code, shaking cache misses. For example, some
7 operators should have iteration first over columns, then rows. ome operators.
8 Since those operators are used in many problems, an implementation must be
9 efficient.
10
11*/
12
13#ifndef __BILINEAR_FORMS_INTEGRATORS_IMPL_HPP__
14#define __BILINEAR_FORMS_INTEGRATORS_IMPL_HPP__
15
16namespace MoFEM {
17
18using CacheMatsTypeType =
19 std::map<std::pair<int, int>, boost::shared_ptr<MatrixDouble>>;
20
21template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
22 typename OpBase>
23struct OpGradGradImpl {};
24
25template <int SPACE_DIM, typename OpBase>
26struct OpGradGradImpl<1, 1, SPACE_DIM, GAUSS, OpBase> : public OpBase {
27 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
28 OpGradGradImpl(const std::string row_field_name,
29 const std::string col_field_name,
30 ScalarFun beta = scalar_fun_one,
31 boost::shared_ptr<Range> ents_ptr = nullptr)
32 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
33 betaCoeff(beta) {
34 if (row_field_name == col_field_name)
35 this->sYmm = true;
36 }
37
38protected:
39 ScalarFun betaCoeff;
40 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
41 EntitiesFieldData::EntData &col_data);
42};
43
44template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
45struct OpGradGradImpl<1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase> : public OpBase {
46 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
47 OpGradGradImpl(const std::string row_field_name,
48 const std::string col_field_name,
49 ScalarFun beta = scalar_fun_one,
50 boost::shared_ptr<Range> ents_ptr = nullptr)
51 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
52 betaCoeff(beta) {
53 if (row_field_name == col_field_name)
54 this->sYmm = true;
55 }
56
57protected:
58 ScalarFun betaCoeff;
59 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
60 EntitiesFieldData::EntData &col_data);
61};
62
63template <int BASE_DIM, int FIELD_DIM, IntegrationType I, typename OpBase>
64struct OpMassImpl {};
65
66template <typename OpBase>
67struct OpMassImpl<1, 1, GAUSS, OpBase> : public OpBase {
68
69 OpMassImpl(
70 const std::string row_field_name, const std::string col_field_name,
71 ScalarFun beta = [](double, double, double) constexpr { return 1; },
72 boost::shared_ptr<Range> ents_ptr = nullptr,
73 boost::shared_ptr<MatrixDouble> cache_mat = nullptr)
74 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
75 betaCoeff(beta) {
76 if (row_field_name == col_field_name)
77 this->sYmm = true;
78 }
79
80protected:
81 using OpBase::OpBase;
82
83 ScalarFun betaCoeff = [](double, double, double) constexpr { return 1; };
84 MoFEMErrorCode integrateImpl(EntitiesFieldData::EntData &row_data,
85 EntitiesFieldData::EntData &col_data,
86 double vol);
87 inline MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
88 EntitiesFieldData::EntData &col_data) {
89 return integrateImpl(row_data, col_data, OpBase::getMeasure());
90 }
91};
92
93template <int FIELD_DIM, typename OpBase>
94struct OpMassImpl<1, FIELD_DIM, GAUSS, OpBase>
95 : public OpMassImpl<1, 1, GAUSS, OpBase> {
96 using OpMassImpl<1, 1, GAUSS, OpBase>::OpMassImpl;
97
98protected:
99 using OpMassImpl<1, 1, GAUSS, OpBase>::betaCoeff;
100 MoFEMErrorCode integrateImpl(EntitiesFieldData::EntData &row_data,
101 EntitiesFieldData::EntData &col_data,
102 double vol);
103 inline MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
104 EntitiesFieldData::EntData &col_data) {
105 return integrateImpl(row_data, col_data, OpBase::getMeasure());
106 }
110};
111
112template <int FIELD_DIM, typename OpBase>
113struct OpMassImpl<3, FIELD_DIM, GAUSS, OpBase> : public OpBase {
114
115 OpMassImpl(const std::string row_field_name, const std::string col_field_name,
116 ScalarFun beta = scalar_fun_one,
117 boost::shared_ptr<Range> ents_ptr = nullptr)
118 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
119 betaCoeff(beta) {
120 if (row_field_name == col_field_name)
121 this->sYmm = true;
122 }
123
124protected:
125 ScalarFun betaCoeff;
126 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
127 EntitiesFieldData::EntData &col_data);
128};
129
130template <typename OpBase>
131struct OpMassImpl<3, 4, GAUSS, OpBase> : public OpBase {
132 OpMassImpl(const std::string row_field_name, const std::string col_field_name,
133 ScalarFun beta = scalar_fun_one,
134 boost::shared_ptr<Range> ents_ptr = nullptr)
135 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
136 betaCoeff(beta) {
137 if (row_field_name == col_field_name)
138 this->sYmm = true;
139 }
140
141protected:
142 ScalarFun betaCoeff;
143 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
144 EntitiesFieldData::EntData &col_data);
145};
146
147template <typename OpBase>
148struct OpMassImpl<3, 9, GAUSS, OpBase> : public OpBase {
149 OpMassImpl(const std::string row_field_name, const std::string col_field_name,
150 ScalarFun beta = scalar_fun_one,
151 boost::shared_ptr<Range> ents_ptr = nullptr)
152 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
153 betaCoeff(beta) {
154 if (row_field_name == col_field_name)
155 this->sYmm = true;
156 }
157
158protected:
159 ScalarFun betaCoeff;
160 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
161 EntitiesFieldData::EntData &col_data);
162};
163
164template <int BASE_DIM, int FIELD_DIM, IntegrationType I, typename OpBase>
165struct OpMassCacheImpl;
166
167template <int FIELD_DIM, IntegrationType I, typename OpBase>
168struct OpMassCacheImpl<1, FIELD_DIM, I, OpBase>
169 : public OpMassImpl<1, FIELD_DIM, I, OpBase> {
170 OpMassCacheImpl(CacheMatsTypeType cache, const std::string row_field_name,
171 const std::string col_field_name, const double beta,
172 boost::shared_ptr<Range> ents_ptr = nullptr)
173 : OpMassImpl<1, FIELD_DIM, I, OpBase>(
174 row_field_name, col_field_name,
175 [](double, double, double) constexpr { return 1; }, ents_ptr),
176 cacheLocMats(cache), scalarBeta(beta) {}
177
178 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
179 EntitiesFieldData::EntData &col_data);
180
181protected:
182 double scalarBeta;
183 CacheMatsTypeType cacheLocMats;
184};
185
186template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
187 typename OpBase>
188struct OpGradSymTensorGradImpl {};
189
190template <int SPACE_DIM, int S, typename OpBase>
191struct OpGradSymTensorGradImpl<1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase>
192 : public OpBase {
193 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
194 OpGradSymTensorGradImpl(const std::string row_field_name,
195 const std::string col_field_name,
196 boost::shared_ptr<MatrixDouble> mat_D,
197 boost::shared_ptr<Range> ents_ptr = nullptr,
198 ScalarFun beta = scalar_fun_one)
199 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
200 matD(mat_D), betaCoeff(beta) {
201 if (row_field_name == col_field_name)
202 this->sYmm = true;
203 }
204
205protected:
206 boost::shared_ptr<MatrixDouble> matD;
207 ScalarFun betaCoeff;
208 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
209 EntitiesFieldData::EntData &col_data);
210};
211
212template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
213 typename OpBase>
214struct OpGradTensorGradImpl {};
215
216template <int SPACE_DIM, int S, typename OpBase>
217struct OpGradTensorGradImpl<1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase>
218 : public OpBase {
219 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
220 OpGradTensorGradImpl(const std::string row_field_name,
221 const std::string col_field_name,
222 boost::shared_ptr<MatrixDouble> mat_D,
223 ScalarFun beta = scalar_fun_one)
224 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL), matD(mat_D),
225 betaCoeff(beta) {}
226
227protected:
228 boost::shared_ptr<MatrixDouble> matD;
229 ScalarFun betaCoeff;
230 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
231 EntitiesFieldData::EntData &col_data);
232};
233
234template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
235 typename OpBase>
236struct OpGradGradSymTensorGradGradImpl {};
237
238template <int FIELD_DIM, int SPACE_DIM, int S, typename OpBase>
239struct OpGradGradSymTensorGradGradImpl<1, FIELD_DIM, SPACE_DIM, S, GAUSS,
240 OpBase> : public OpBase {
241 OpGradGradSymTensorGradGradImpl(const std::string row_field_name,
242 const std::string col_field_name,
243 boost::shared_ptr<MatrixDouble> mat_D,
244 boost::shared_ptr<Range> ents_ptr = nullptr)
245 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
246 matD(mat_D) {
247 if (row_field_name == col_field_name)
248 this->sYmm = true;
249 }
250
251protected:
252 boost::shared_ptr<MatrixDouble> matD;
253 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
254 EntitiesFieldData::EntData &col_data);
255};
256
257template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, int S, IntegrationType I,
258 typename OpBase>
259struct OpGradGradTensorGradGradImpl {};
260
261template <int FIELD_DIM, int SPACE_DIM, int S, typename OpBase>
262struct OpGradGradTensorGradGradImpl<1, FIELD_DIM, SPACE_DIM, S, GAUSS, OpBase>
263 : public OpBase {
264 OpGradGradTensorGradGradImpl(const std::string row_field_name,
265 const std::string col_field_name,
266 boost::shared_ptr<MatrixDouble> mat_D,
267 boost::shared_ptr<Range> ents_ptr = nullptr)
268 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
269 matD(mat_D) {
270 if (row_field_name == col_field_name)
271 this->sYmm = true;
272 }
273
274protected:
275 boost::shared_ptr<MatrixDouble> matD;
276 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
277 EntitiesFieldData::EntData &col_data);
278};
279
280template <int SPACE_DIM, IntegrationType I, typename OpBase>
281struct OpMixDivTimesScalarImpl {};
282
283template <int SPACE_DIM, typename OpBase>
284struct OpMixDivTimesScalarImpl<SPACE_DIM, GAUSS, OpBase> : public OpBase {
285 OpMixDivTimesScalarImpl(const std::string row_field_name,
286 const std::string col_field_name,
287 ConstantFun alpha_fun,
288 const bool assemble_transpose = false,
289 const bool only_transpose = false)
290 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
291 alphaConstant(alpha_fun) {
292 this->assembleTranspose = assemble_transpose;
293 this->onlyTranspose = only_transpose;
294 }
295
296protected:
297 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
298 ConstantFun alphaConstant;
299 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
300 EntitiesFieldData::EntData &col_data);
301};
302
303template <int SPACE_DIM, IntegrationType I, typename OpBase,
304 CoordinateTypes CoordSys>
305struct OpMixDivTimesVecImpl {};
306
307template <int SPACE_DIM, typename OpBase, CoordinateTypes CoordSys>
308struct OpMixDivTimesVecImpl<SPACE_DIM, GAUSS, OpBase, CoordSys>
309 : public OpBase {
310
311 OpMixDivTimesVecImpl(const std::string row_field_name,
312 const std::string col_field_name, ConstantFun alpha_fun,
313 const bool assemble_transpose,
314 const bool only_transpose = false)
315 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
316 alphaConstant(alpha_fun) {
317 this->assembleTranspose = assemble_transpose;
318 this->onlyTranspose = only_transpose;
319 }
320
321 OpMixDivTimesVecImpl(const std::string row_field_name,
322 const std::string col_field_name, ConstantFun alpha_fun,
323 ScalarFun beta_fun, const bool assemble_transpose,
324 const bool only_transpose = false)
325 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
326 alphaConstant(alpha_fun), betaCoeff(beta_fun) {
327 this->assembleTranspose = assemble_transpose;
328 this->onlyTranspose = only_transpose;
329 }
330
331protected:
332 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
333
334 ConstantFun alphaConstant = []() constexpr { return 1; };
335 ScalarFun betaCoeff = [](double, double, double) constexpr { return 1; };
336
337 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
338 EntitiesFieldData::EntData &col_data);
339};
340
341template <int SPACE_DIM, IntegrationType I, typename OpBase,
342 CoordinateTypes COORDINATE_SYSTEM>
343struct OpMixScalarTimesDivImpl {};
344
345template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
346struct OpMixScalarTimesDivImpl<SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM>
347 : public OpBase {
348 OpMixScalarTimesDivImpl(
349 const std::string row_field_name, const std::string col_field_name,
350 ScalarFun alpha_fun = [](double, double, double) constexpr { return 1; },
351 const bool assemble_transpose = false, const bool only_transpose = false)
352 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
353 alphaConstant(alpha_fun) {
354 this->assembleTranspose = assemble_transpose;
355 this->onlyTranspose = only_transpose;
356 this->sYmm = false;
357 }
358
359protected:
360 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
361 ScalarFun alphaConstant;
362 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
363 EntitiesFieldData::EntData &col_data);
364};
365
366template <int BASE_DIM, int FIELD_DIM, int SPACE_DIM, IntegrationType I,
367 typename OpBase>
368struct OpMixVectorTimesGradImpl;
369
370template <int SPACE_DIM, typename OpBase>
371struct OpMixVectorTimesGradImpl<3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase>
372 : public OpBase {
373 OpMixVectorTimesGradImpl(const std::string row_field_name,
374 const std::string col_field_name,
375 ConstantFun alpha_fun,
376 const bool assemble_transpose = false,
377 const bool only_transpose = false)
378 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
379 alphaConstant(alpha_fun) {
380 this->assembleTranspose = assemble_transpose;
381 this->onlyTranspose = only_transpose;
382 }
383
384protected:
385 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
386 ConstantFun alphaConstant;
387 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
388 EntitiesFieldData::EntData &col_data);
389};
390
391template <int SPACE_DIM, typename OpBase>
392struct OpMixVectorTimesGradImpl<1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase>
393 : public OpBase {
394 OpMixVectorTimesGradImpl(const std::string row_field_name,
395 const std::string col_field_name,
396 ConstantFun alpha_fun,
397 const bool assemble_transpose = false,
398 const bool only_transpose = false)
399 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
400 alphaConstant(alpha_fun) {
401 this->assembleTranspose = assemble_transpose;
402 this->onlyTranspose = only_transpose;
403 }
404
405protected:
406 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
407 ConstantFun alphaConstant;
408 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
409 EntitiesFieldData::EntData &col_data);
410};
411
412template <int SPACE_DIM, IntegrationType I, typename OpBase>
413struct OpMixTensorTimesGradImpl {};
414
415template <int SPACE_DIM, typename OpBase>
416struct OpMixTensorTimesGradImpl<SPACE_DIM, GAUSS, OpBase> : public OpBase {
417
418 OpMixTensorTimesGradImpl(const std::string row_field_name,
419 const std::string col_field_name,
420 ConstantFun alpha_fun, const bool assemble_transpose,
421 const bool only_transpose = false)
422 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
423 alphaConstant(alpha_fun) {
424 this->assembleTranspose = assemble_transpose;
425 this->onlyTranspose = only_transpose;
426 }
427
428 OpMixTensorTimesGradImpl(const std::string row_field_name,
429 const std::string col_field_name,
430 ConstantFun alpha_fun, ScalarFun beta_coeff,
431 const bool assemble_transpose,
432 const bool only_transpose = false)
433 : OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
434 alphaConstant(alpha_fun), betaCoeff(beta_coeff) {
435 this->assembleTranspose = assemble_transpose;
436 this->onlyTranspose = only_transpose;
437 }
438
439protected:
440 FTensor::Index<'i', SPACE_DIM> i; ///< summit Index
441 FTensor::Index<'j', SPACE_DIM> j; ///< summit Index
442 ConstantFun alphaConstant = []() constexpr { return 1; };
443 ScalarFun betaCoeff = [](double, double, double) constexpr { return 1; };
444 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
445 EntitiesFieldData::EntData &col_data);
446};
447
448template <int SPACE_DIM, typename OpBase>
449MoFEMErrorCode OpGradGradImpl<1, 1, SPACE_DIM, GAUSS, OpBase>::iNtegrate(
450 EntitiesFieldData::EntData &row_data,
451 EntitiesFieldData::EntData &col_data) {
452 MoFEMFunctionBegin;
453
454 // get element volume
455 const double vol = OpBase::getMeasure();
456 // get integration weights
457 auto t_w = OpBase::getFTensor0IntegrationWeight();
458 // get base function gradient on rows
459 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
460 // get coordinate at integration points
461 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
462
463 // loop over integration points
464 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
465 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
466 // take into account Jacobian
467 const double alpha = t_w * beta;
468 // loop over ros base functions
469 int rr = 0;
470 for (; rr != OpBase::nbRows; rr++) {
471 // get column base functions gradient at gauss point gg
472 auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
473 // loop over columns
474 for (int cc = 0; cc != OpBase::nbCols; cc++) {
475 // calculate element of local matrix
476 OpBase::locMat(rr, cc) += alpha * (t_row_grad(i) * t_col_grad(i));
477 ++t_col_grad; // move to another gradient of base function
478 // on column
479 }
480 ++t_row_grad; // move to another element of gradient of base
481 // function on row
482 }
483 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
484 ++t_row_grad;
485
486 ++t_coords;
487 ++t_w; // move to another integration weight
488 }
489 MoFEMFunctionReturn(0);
490}
491
492template <int FIELD_DIM, int SPACE_DIM, typename OpBase>
493MoFEMErrorCode
494OpGradGradImpl<1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase>::iNtegrate(
495 EntitiesFieldData::EntData &row_data,
496 EntitiesFieldData::EntData &col_data) {
497 MoFEMFunctionBegin;
498
499 // get element volume
500 const double vol = OpBase::getMeasure();
501 // get integration weights
502 auto t_w = OpBase::getFTensor0IntegrationWeight();
503 // get base function gradient on rows
504 auto t_row_grad = row_data.getFTensor1DiffN<SPACE_DIM>();
505 // get coordinate at integration points
506 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
507
508 FTensor::Index<'i', SPACE_DIM> i;
509 FTensor::Index<'j', FIELD_DIM> j;
510
511 auto get_t_vec = [&](const int rr) {
512 std::array<double *, FIELD_DIM> ptrs;
513 for (auto i = 0; i != FIELD_DIM; ++i)
514 ptrs[i] = &OpBase::locMat(rr + i, i);
515 return FTensor::Tensor1<FTensor::PackPtr<double *, FIELD_DIM>, FIELD_DIM>(
516 ptrs);
517 };
518
519 // loop over integration points
520 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
521 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
522 // take into account Jacobian
523 const double alpha = t_w * beta;
524 // loop over ros base functions
525 int rr = 0;
526 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
527 // get diag vec
528 auto t_vec = get_t_vec(rr * FIELD_DIM);
529 // get column base functions gradient at gauss point gg
530 auto t_col_grad = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
531 // loop over columns
532 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
533 // calculate element of local matrix
534 t_vec(j) += alpha * (t_row_grad(i) * t_col_grad(i));
535 ++t_col_grad; // move to another gradient of base function
536 // on column
537 ++t_vec;
538 }
539 ++t_row_grad; // move to another element of gradient of base
540 // function on row
541 }
542 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
543 ++t_row_grad;
544
545 ++t_coords;
546 ++t_w; // move to another integration weight
547 }
548 MoFEMFunctionReturn(0);
549}
550
551template <typename OpBase>
552MoFEMErrorCode OpMassImpl<1, 1, GAUSS, OpBase>::integrateImpl(
553 EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data,
554 double vol) {
555 MoFEMFunctionBegin;
556
557#ifndef NDEBUG
558 auto log_error = [&]() {
559 MOFEM_LOG("SELF", Sev::error) << "Row side " << OpBase::rowSide << " "
560 << CN::EntityTypeName(OpBase::rowType);
561 MOFEM_LOG("SELF", Sev::error) << "Col side " << OpBase::colSide << " "
562 << CN::EntityTypeName(OpBase::colType);
563 };
564
565 if (row_data.getN().size2() < OpBase::nbRows) {
566 log_error();
567 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
568 "Wrong number of base functions on rows %zu < %d",
569 row_data.getN().size2(), OpBase::nbRows);
570 }
571 if (col_data.getN().size2() < OpBase::nbCols) {
572 log_error();
573 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
574 "Wrong number of base functions on cols %zu < %d",
575 col_data.getN().size2(), OpBase::nbCols);
576 }
577 if (row_data.getN().size1() != OpBase::nbIntegrationPts) {
578 log_error();
579 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
580 "Wrong number of integration points on rows %zu != %d",
581 row_data.getN().size1(), OpBase::nbIntegrationPts);
582 }
583 if (col_data.getN().size1() != OpBase::nbIntegrationPts) {
584 log_error();
585 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
586 "Wrong number of integration points on cols %zu != %d",
587 col_data.getN().size1(), OpBase::nbIntegrationPts);
588 }
589#endif
590
591 // get integration weights
592 auto t_w = OpBase::getFTensor0IntegrationWeight();
593 // get base function gradient on rows
594 auto t_row_base = row_data.getFTensor0N();
595 // get coordinate at integration points
596 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
597 // loop over integration points
598 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
599 const double beta = betaCoeff(t_coords(0), t_coords(1), t_coords(2));
600 // take into account Jacobian
601 const double alpha = t_w * beta;
602 // loop over rows base functions
603 auto a_mat_ptr = &*OpBase::locMat.data().begin();
604 int rr = 0;
605 for (; rr != OpBase::nbRows; rr++) {
606 // get column base functions gradient at gauss point gg
607 auto t_col_base = col_data.getFTensor0N(gg, 0);
608 // loop over columns
609 for (int cc = 0; cc != OpBase::nbCols; cc++) {
610 // calculate element of local matrix
611 *a_mat_ptr += alpha * (t_row_base * t_col_base);
612 ++t_col_base;
613 ++a_mat_ptr;
614 }
615 ++t_row_base;
616 }
617 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
618 ++t_row_base;
619 ++t_coords;
620 ++t_w; // move to another integration weight
621 }
622
623 OpBase::locMat *= vol;
624
625 MoFEMFunctionReturn(0);
626};
627
628template <int FIELD_DIM, typename OpBase>
629MoFEMErrorCode OpMassImpl<1, FIELD_DIM, GAUSS, OpBase>::integrateImpl(
630 EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data,
631 double vol) {
632 MoFEMFunctionBegin;
633
634 auto integrate = [&](auto &mat) {
635 MoFEMFunctionBeginHot;
636 // get integration weights
637 auto t_w = OpBase::getFTensor0IntegrationWeight();
638 // get base function gradient on rows
639 auto t_row_base = row_data.getFTensor0N();
640 // get coordinate at integration points
641 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
642
643 FTensor::Index<'i', FIELD_DIM> i;
644 auto get_t_vec = [&](const int rr) {
645 std::array<double *, FIELD_DIM> ptrs;
646 for (auto i = 0; i != FIELD_DIM; ++i)
647 ptrs[i] = &mat(rr + i, i);
648 return FTensor::Tensor1<FTensor::PackPtr<double *, FIELD_DIM>, FIELD_DIM>(
649 ptrs);
650 };
651
652 // loop over integration points
653 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
654 const double beta = betaCoeff(t_coords(0), t_coords(1), t_coords(2));
655 // take into account Jacobian
656 const double alpha = t_w * beta;
657 // loop over rows base functions
658 int rr = 0;
659 for (; rr != OpBase::nbRows / FIELD_DIM; rr++) {
660 // get column base functions gradient at gauss point gg
661 auto t_col_base = col_data.getFTensor0N(gg, 0);
662 // get mat vec
663 auto t_vec = get_t_vec(FIELD_DIM * rr);
664 // loop over columns
665 for (int cc = 0; cc != OpBase::nbCols / FIELD_DIM; cc++) {
666 // calculate element of local matrix
667 t_vec(i) += alpha * (t_row_base * t_col_base);
668 ++t_col_base;
669 ++t_vec;
670 }
671 ++t_row_base;
672 }
673 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
674 ++t_row_base;
675 ++t_coords;
676 ++t_w; // move to another integration weight
677 }
678 MoFEMFunctionReturnHot(0);
679 };
680
681 CHKERR integrate(OpBase::locMat);
682 OpBase::locMat *= vol;
683
684 MoFEMFunctionReturn(0);
685};
686
687template <int FIELD_DIM, typename OpBase>
688MoFEMErrorCode OpMassImpl<3, FIELD_DIM, GAUSS, OpBase>::iNtegrate(
689 EntitiesFieldData::EntData &row_data,
690 EntitiesFieldData::EntData &col_data) {
691 FTensor::Index<'i', FIELD_DIM> i;
692 MoFEMFunctionBegin;
693 size_t nb_base_functions = row_data.getN().size2() / 3;
694 // // get element volume
695 const double vol = OpBase::getMeasure();
696 // get integration weights
697 auto t_w = OpBase::getFTensor0IntegrationWeight();
698 // get base function gradient on rows
699 auto t_row_base = row_data.getFTensor1N<3>();
700 // get coordinate at integration points
701 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
702 // loop over integration points
703 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
704 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
705 // take into account Jacobian
706 const double alpha = t_w * beta;
707 // loop over rows base functions
708 auto a_mat_ptr = &*OpBase::locMat.data().begin();
709 int rr = 0;
710 for (; rr != OpBase::nbRows; rr++) {
711 // get column base functions gradient at gauss point gg
712 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
713 // loop over columns
714 for (int cc = 0; cc != OpBase::nbCols; cc++) {
715 // calculate element of local matrix
716 (*a_mat_ptr) += alpha * (t_row_base(i) * t_col_base(i));
717 ++t_col_base;
718 ++a_mat_ptr;
719 }
720 ++t_row_base;
721 }
722 for (; rr < nb_base_functions; ++rr)
723 ++t_row_base;
724 ++t_coords;
725 ++t_w; // move to another integration weight
726 }
727 MoFEMFunctionReturn(0);
728};
729
730template <typename OpBase>
731MoFEMErrorCode OpMassImpl<3, 4, GAUSS, OpBase>::iNtegrate(
732 EntitiesFieldData::EntData &row_data,
733 EntitiesFieldData::EntData &col_data) {
734 MoFEMFunctionBegin;
735 FTensor::Index<'i', 2> I;
736 FTensor::Index<'k', 3> k;
737 auto get_t_vec = [&](const int rr) {
738 return FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{
739 &OpBase::locMat(rr + 0, 0), &OpBase::locMat(rr + 1, 1)};
740 };
741 size_t nb_base_functions = row_data.getN().size2() / 3;
742 // // get element volume
743 const double vol = OpBase::getMeasure();
744 // get integration weights
745 auto t_w = OpBase::getFTensor0IntegrationWeight();
746 // get base function gradient on rows
747 auto t_row_base = row_data.getFTensor1N<3>();
748 // get coordinate at integration points
749 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
750 // loop over integration points
751 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
752 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
753 // take into account Jacobian
754 const double alpha = t_w * beta;
755 // loop over rows base functions
756 int rr = 0;
757 for (; rr != OpBase::nbRows / 2; rr++) {
758 // get column base functions gradient at gauss point gg
759 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
760 auto t_vec = get_t_vec(2 * rr);
761 // loop over columns
762 for (int cc = 0; cc != OpBase::nbCols / 2; cc++) {
763 // calculate element of local matrix
764 t_vec(I) += alpha * (t_row_base(k) * t_col_base(k));
765 ++t_col_base;
766 ++t_vec;
767 }
768 ++t_row_base;
769 }
770 for (; rr < nb_base_functions; ++rr)
771 ++t_row_base;
772 ++t_coords;
773 ++t_w; // move to another integration weight
774 }
775 MoFEMFunctionReturn(0);
776}
777
778template <typename OpBase>
779MoFEMErrorCode OpMassImpl<3, 9, GAUSS, OpBase>::iNtegrate(
780 EntitiesFieldData::EntData &row_data,
781 EntitiesFieldData::EntData &col_data) {
782 MoFEMFunctionBegin;
783 FTensor::Index<'i', 3> i;
784 FTensor::Index<'k', 3> k;
785 auto get_t_vec = [&](const int rr) {
786 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>{
787 &OpBase::locMat(rr + 0, 0), &OpBase::locMat(rr + 1, 1),
788 &OpBase::locMat(rr + 2, 2)};
789 };
790 size_t nb_base_functions = row_data.getN().size2() / 3;
791 // // get element volume
792 const double vol = OpBase::getMeasure();
793 // get integration weights
794 auto t_w = OpBase::getFTensor0IntegrationWeight();
795 // get base function gradient on rows
796 auto t_row_base = row_data.getFTensor1N<3>();
797 // get coordinate at integration points
798 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
799 // loop over integration points
800 for (int gg = 0; gg != OpBase::nbIntegrationPts; gg++) {
801 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
802 // take into account Jacobian
803 const double alpha = t_w * beta;
804 // loop over rows base functions
805 int rr = 0;
806 for (; rr != OpBase::nbRows / 3; rr++) {
807 // get column base functions gradient at gauss point gg
808 auto t_col_base = col_data.getFTensor1N<3>(gg, 0);
809 auto t_vec = get_t_vec(3 * rr);
810 // loop over columns
811 for (int cc = 0; cc != OpBase::nbCols / 3; cc++) {
812 // calculate element of local matrix
813 t_vec(i) += alpha * (t_row_base(k) * t_col_base(k));
814 ++t_col_base;
815 ++t_vec;
816 }
817 ++t_row_base;
818 }
819 for (; rr < nb_base_functions; ++rr)
820 ++t_row_base;
821 ++t_coords;
822 ++t_w; // move to another integration weight
823 }
824 MoFEMFunctionReturn(0);
825};
826
827template <int FIELD_DIM, IntegrationType I, typename OpBase>
828MoFEMErrorCode OpMassCacheImpl<1, FIELD_DIM, I, OpBase>::iNtegrate(
829 EntitiesFieldData::EntData &row_data,
830 EntitiesFieldData::EntData &col_data) {
831 MoFEMFunctionBegin;
832
833 const auto vol = this->getMeasure();
834 const auto row_type = this->rowType;
835 const auto col_type = this->colType;
836 auto &loc_mat = this->locMat;
837
838 auto p = std::make_pair(row_type, col_type);
839
840 if (cacheLocMats[p]) {
841 if (cacheLocMats[p]->size1() != loc_mat.size1() &&
842 cacheLocMats[p]->size2() != loc_mat.size2()) {
843 cacheLocMats[p]->resize(loc_mat.size1(), loc_mat.size2());
844 CHKERR this->integrateImpl(row_data, col_data, 1);
845 *(cacheLocMats[p]) = loc_mat;
846 } else {
847 loc_mat = *(cacheLocMats[p]);
848 }
849 loc_mat *= scalarBeta * this->getMeasure();
850 } else {
851 CHKERR this->integrateImpl(row_data, col_data,
852 scalarBeta * this->getMeasure());
853 }
854 MoFEMFunctionReturn(0);
855}
856
857template <int SPACE_DIM, int S, typename OpBase>
858MoFEMErrorCode
859OpGradSymTensorGradImpl<1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase>::iNtegrate(
860 EntitiesFieldData::EntData &row_data,
861 EntitiesFieldData::EntData &col_data) {
862 MoFEMFunctionBegin;
863
864 const size_t nb_row_dofs = row_data.getIndices().size();
865 const size_t nb_col_dofs = col_data.getIndices().size();
866
867 if (nb_row_dofs && nb_col_dofs) {
868
869 const bool diag = (row_data.getFieldEntities()[0]->getLocalUniqueId() ==
870 col_data.getFieldEntities()[0]->getLocalUniqueId());
871
872 FTensor::Index<'i', SPACE_DIM> i;
873 FTensor::Index<'j', SPACE_DIM> j;
874 FTensor::Index<'k', SPACE_DIM> k;
875 FTensor::Index<'l', SPACE_DIM> l;
876
877 // get element volume
878 double vol = OpBase::getMeasure();
879
880 // get intergrayion weights
881 auto t_w = OpBase::getFTensor0IntegrationWeight();
882
883 // get derivatives of base functions on rows
884 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
885
886 // Elastic stiffness tensor (4th rank tensor with minor and major
887 // symmetry)
888 using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
889 auto get_D_at_pts =
890 MatrixSizeHelper<GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, S, DL>,
891 DL>::get(*matD, OpBase::nbIntegrationPts);
892 auto t_D_at_pts = get_D_at_pts();
893
894 // get coordinate at integration points
895 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
896
897 // iterate over integration points
898 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
899
900 // calculate scalar weight times element volume
901 double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
902
903 // iterate over row base functions
904 int rr = 0;
905 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
906
907 // get sub matrix for the row
908 auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
909
910 FTensor::Christof<double, SPACE_DIM, SPACE_DIM> t_rowD;
911 // I mix up the indices here so that it behaves like a
912 // Dg. That way I don't have to have a separate wrapper
913 // class Christof_Expr, which simplifies things.
914 t_rowD(l, j, k) = t_D_at_pts(i, j, k, l) * (a * t_row_diff_base(i));
915
916 // get derivatives of base functions for columns
917 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
918
919 // iterate column base functions
920 const auto nb_cols = (diag) ? rr : OpBase::nbCols / SPACE_DIM - 1;
921 for (int cc = 0; cc <= nb_cols; ++cc) {
922
923 // integrate block local stiffens matrix
924 t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
925
926 // move to next column base function
927 ++t_col_diff_base;
928
929 // move to next block of local stiffens matrix
930 ++t_m;
931 }
932
933 // move to next row base function
934 ++t_row_diff_base;
935 }
936
937 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
938 ++t_row_diff_base;
939
940 // move to next integration weight
941 ++t_w;
942 ++t_D_at_pts;
943 ++t_coords;
944 }
945
946 // Copy symmetry
947 if (diag) {
948 for (int rr = 0; rr != OpBase::nbRows / SPACE_DIM - 1; ++rr) {
949 auto t_m_rr = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
950 this->locMat, SPACE_DIM * rr, SPACE_DIM * (rr + 1));
951 auto t_m_cc = getFTensor2FromArray<SPACE_DIM, SPACE_DIM>(
952 this->locMat, SPACE_DIM * (rr + 1), SPACE_DIM * rr,
953 SPACE_DIM * OpBase::nbCols);
954 for (int cc = rr + 1; cc < OpBase::nbCols / SPACE_DIM; ++cc) {
955 t_m_rr(i, j) = t_m_cc(j, i);
956 ++t_m_rr;
957 ++t_m_cc;
958 }
959 }
960 }
961 }
962
963 MoFEMFunctionReturn(0);
964}
965
966template <int FIELD_DIM, int SPACE_DIM, int S, typename OpBase>
967MoFEMErrorCode
968OpGradGradSymTensorGradGradImpl<1, FIELD_DIM, SPACE_DIM, S, GAUSS, OpBase>::
969 iNtegrate(EntitiesFieldData::EntData &row_data,
970 EntitiesFieldData::EntData &col_data) {
971 MoFEMFunctionBegin;
972
973 if (OpBase::nbRows && OpBase::nbCols) {
974
975 FTensor::Index<'i', SPACE_DIM> i;
976 FTensor::Index<'j', SPACE_DIM> j;
977 FTensor::Index<'k', SPACE_DIM> k;
978 FTensor::Index<'l', SPACE_DIM> l;
979 FTensor::Index<'I', FIELD_DIM> I;
980 FTensor::Index<'J', FIELD_DIM> J;
981
982 auto &row_hessian = row_data.getN(BaseDerivatives::SecondDerivative);
983 auto &col_hessian = col_data.getN(BaseDerivatives::SecondDerivative);
984
985#ifndef NDEBUG
986 if (row_hessian.size1() != OpBase::nbIntegrationPts) {
987 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
988 "Wrong number of integration pts (%zu != %d)", row_hessian.size1(),
989 OpBase::nbIntegrationPts);
990 }
991 if (row_hessian.size2() !=
992 OpBase::nbRowBaseFunctions * SPACE_DIM * SPACE_DIM) {
993 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
994 "Wrong number of base functions (%zu != %d)",
995 row_hessian.size2() / (SPACE_DIM * SPACE_DIM),
996 OpBase::nbRowBaseFunctions);
997 }
998 if (row_hessian.size2() < OpBase::nbRows * SPACE_DIM * SPACE_DIM) {
999 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1000 "Wrong number of base functions (%zu < %d)", row_hessian.size2(),
1001 OpBase::nbRows * SPACE_DIM * SPACE_DIM);
1002 }
1003 if (col_hessian.size1() != OpBase::nbIntegrationPts) {
1004 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1005 "Wrong number of integration pts (%zu != %d)", col_hessian.size1(),
1006 OpBase::nbIntegrationPts);
1007 }
1008 if (col_hessian.size2() < OpBase::nbCols * SPACE_DIM * SPACE_DIM) {
1009 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1010 "Wrong number of base functions (%zu < %d)", col_hessian.size2(),
1011 OpBase::nbCols * SPACE_DIM * SPACE_DIM);
1012 }
1013#endif
1014
1015 // get element volume
1016 double vol = OpBase::getMeasure();
1017
1018 // get intergrayion weights
1019 auto t_w = OpBase::getFTensor0IntegrationWeight();
1020
1021 auto t_row_diff2 = getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(
1022 &*row_hessian.data().begin());
1023
1024 // Elastic stiffness tensor (4th rank tensor with minor and major
1025 // symmetry)
1026 using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
1027 auto get_D_at_pts =
1028 MatrixSizeHelper<GetFTensor4DdgFromMatType<SPACE_DIM, SPACE_DIM, S, DL>,
1029 DL>::get(*matD, OpBase::nbIntegrationPts);
1030 auto t_D_at_pts = get_D_at_pts();
1031 auto t_kd = FTensor::Kronecker_Delta<double>();
1032
1033 // iterate over integration points
1034 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1035
1036 // calculate scalar weight times element volume
1037 double a = t_w * vol;
1038
1039 // iterate over row base functions
1040 int rr = 0;
1041 for (; rr != OpBase::nbRows; ++rr) {
1042
1043 auto t_mat = getFTensor2FromPtr<FIELD_DIM, FIELD_DIM, FIELD_DIM>(
1044 &OpBase::locMat(rr, 0));
1045
1046 FTensor::Tensor2_symmetric<double, SPACE_DIM> t_rowD;
1047 t_rowD(k, l) = t_D_at_pts(i, j, k, l) * (a * t_row_diff2(i, j));
1048
1049 // get derivatives of base functions for columns
1050 auto t_col_diff2 =
1051 getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(&col_hessian(gg, 0));
1052
1053 // iterate column base functions
1054 for (int cc = 0; cc != OpBase::nbCols; ++cc) {
1055
1056 // integrate block local stiffens matrix
1057 t_mat(I, J) += (t_rowD(i, j) * t_col_diff2(i, j)) * t_kd(I, J);
1058
1059 // move to next column base function
1060 ++t_col_diff2;
1061
1062 // move to next block of local stiffens matrix
1063 ++t_mat;
1064 }
1065
1066 // move to next row base function
1067 ++t_row_diff2;
1068 }
1069
1070 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1071 ++t_row_diff2;
1072
1073 // move to next integration weight
1074 ++t_w;
1075 ++t_D_at_pts;
1076 }
1077 }
1078
1079 MoFEMFunctionReturn(0);
1080}
1081
1082template <int FIELD_DIM, int SPACE_DIM, int S, typename OpBase>
1083MoFEMErrorCode
1084OpGradGradTensorGradGradImpl<1, FIELD_DIM, SPACE_DIM, S, GAUSS, OpBase>::
1085 iNtegrate(EntitiesFieldData::EntData &row_data,
1086 EntitiesFieldData::EntData &col_data) {
1087 MoFEMFunctionBegin;
1088
1089 if (OpBase::nbRows && OpBase::nbCols) {
1090
1091 FTensor::Index<'i', SPACE_DIM> i;
1092 FTensor::Index<'j', SPACE_DIM> j;
1093 FTensor::Index<'k', SPACE_DIM> k;
1094 FTensor::Index<'l', SPACE_DIM> l;
1095 FTensor::Index<'I', FIELD_DIM> I;
1096 FTensor::Index<'J', FIELD_DIM> J;
1097
1098 auto &row_hessian = row_data.getN(BaseDerivatives::SecondDerivative);
1099 auto &col_hessian = col_data.getN(BaseDerivatives::SecondDerivative);
1100
1101#ifndef NDEBUG
1102 if (row_hessian.size1() != OpBase::nbIntegrationPts) {
1103 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1104 "Wrong number of integration pts (%zu != %d)", row_hessian.size1(),
1105 OpBase::nbIntegrationPts);
1106 }
1107 if (row_hessian.size2() !=
1108 OpBase::nbRowBaseFunctions * SPACE_DIM * SPACE_DIM) {
1109 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1110 "Wrong number of base functions (%zu != %d)",
1111 row_hessian.size2() / (SPACE_DIM * SPACE_DIM),
1112 OpBase::nbRowBaseFunctions);
1113 }
1114 if (row_hessian.size2() < OpBase::nbRows * SPACE_DIM * SPACE_DIM) {
1115 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1116 "Wrong number of base functions (%zu < %d)", row_hessian.size2(),
1117 OpBase::nbRows * SPACE_DIM * SPACE_DIM);
1118 }
1119 if (col_hessian.size1() != OpBase::nbIntegrationPts) {
1120 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1121 "Wrong number of integration pts (%zu != %d)", col_hessian.size1(),
1122 OpBase::nbIntegrationPts);
1123 }
1124 if (col_hessian.size2() < OpBase::nbCols * SPACE_DIM * SPACE_DIM) {
1125 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1126 "Wrong number of base functions (%zu < %d)", col_hessian.size2(),
1127 OpBase::nbCols * SPACE_DIM * SPACE_DIM);
1128 }
1129#endif
1130
1131 // get element volume
1132 double vol = OpBase::getMeasure();
1133
1134 // get intergrayion weights
1135 auto t_w = OpBase::getFTensor0IntegrationWeight();
1136
1137 auto t_row_diff2 = getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(
1138 &*row_hessian.data().begin());
1139
1140 // Elastic stiffness tensor
1141 using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
1142 auto get_D_at_pts = MatrixSizeHelper<
1143 GetFTensor4FromMatType<SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM, S,
1144 DL>,
1145 DL>::get(*matD, OpBase::nbIntegrationPts);
1146 auto t_D_at_pts = get_D_at_pts();
1147
1148 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
1149
1150 // iterate over integration points
1151 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1152
1153 // calculate scalar weight times element volume
1154 double a = t_w * vol;
1155
1156 // get sub matrix for the row
1157 auto m_ptr = &*OpBase::locMat.data().begin();
1158
1159 // iterate over row base functions
1160 int rr = 0;
1161 for (; rr != OpBase::nbRows; ++rr) {
1162
1163 FTensor::Tensor2_symmetric<double, SPACE_DIM> t_rowD;
1164 t_rowD(k, l) = t_D_at_pts(i, j, k, l) * (a * t_row_diff2(i, j));
1165
1166 // get derivatives of base functions for columns
1167 auto t_col_diff2 =
1168 getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(&col_hessian(gg, 0));
1169
1170 auto t_mat = getFTensor2FromPtr<FIELD_DIM, FIELD_DIM, FIELD_DIM>(
1171 &OpBase::locMat(FIELD_DIM * rr, 0));
1172
1173 // iterate column base functions
1174 for (int cc = 0; cc != OpBase::nbCols; ++cc) {
1175
1176 // integrate block local stiffens matrix
1177 t_mat(I, J) += (t_rowD(i, j) * t_col_diff2(i, j)) * t_kd(I, J);
1178
1179 // move to next column base function
1180 ++t_col_diff2;
1181
1182 // move to next block of local stiffens matrix
1183 ++t_mat;
1184 }
1185
1186 // move to next row base function
1187 ++t_row_diff2;
1188 }
1189
1190 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1191 ++t_row_diff2;
1192
1193 // move to next integration weight
1194 ++t_w;
1195 ++t_D_at_pts;
1196 }
1197 }
1198
1199 MoFEMFunctionReturn(0);
1200}
1201
1202template <int SPACE_DIM, int S, typename OpBase>
1203MoFEMErrorCode
1204OpGradTensorGradImpl<1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase>::iNtegrate(
1205 EntitiesFieldData::EntData &row_data,
1206 EntitiesFieldData::EntData &col_data) {
1207 MoFEMFunctionBegin;
1208
1209 const size_t nb_row_dofs = row_data.getIndices().size();
1210 const size_t nb_col_dofs = col_data.getIndices().size();
1211
1212 if (nb_row_dofs && nb_col_dofs) {
1213
1214 FTensor::Index<'i', SPACE_DIM> i;
1215 FTensor::Index<'j', SPACE_DIM> j;
1216 FTensor::Index<'k', SPACE_DIM> k;
1217 FTensor::Index<'l', SPACE_DIM> l;
1218
1219 // get element volume
1220 double vol = OpBase::getMeasure();
1221
1222 // get intergrayion weights
1223 auto t_w = OpBase::getFTensor0IntegrationWeight();
1224
1225 // get derivatives of base functions on rows
1226 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
1227
1228 // stiffness tensor (4th rank tensor)
1229 using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
1230 auto get_D_at_pts = MatrixSizeHelper<
1231 GetFTensor4FromMatType<SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM, S,
1232 DL>,
1233 DL>::get(*matD, OpBase::nbIntegrationPts);
1234 auto t_D_at_pts = get_D_at_pts();
1235
1236 // get coordinate at integration points
1237 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1238
1239 // iterate over integration points
1240 for (int gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1241
1242 // calculate scalar weight times element volume
1243 double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1244
1245 // iterate over row base functions
1246 int rr = 0;
1247 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1248
1249 // get sub matrix for the row
1250 auto t_m = OpBase::template getLocMat<SPACE_DIM>(SPACE_DIM * rr);
1251
1252 // calculate row
1253 FTensor::Tensor3<double, 3, 3, 3> t_row;
1254 t_row(i, k, l) = t_D_at_pts(i, j, k, l) * (a * t_row_diff_base(j));
1255
1256 // get derivatives of base functions for columns
1257 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1258
1259 // iterate column base functions
1260 for (int cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1261
1262 // integrate block local stiffens matrix
1263 t_m(i, k) += t_row(i, k, l) * t_col_diff_base(l);
1264
1265 // move to next column base function
1266 ++t_col_diff_base;
1267
1268 // move to next block of local stiffens matrix
1269 ++t_m;
1270 }
1271
1272 // move to next row base function
1273 ++t_row_diff_base;
1274 }
1275
1276 for (; rr < OpBase::nbRowBaseFunctions; ++rr)
1277 ++t_row_diff_base;
1278
1279 // move to next integration weight
1280 ++t_w;
1281 ++t_D_at_pts;
1282 ++t_coords;
1283 }
1284 }
1285
1286 MoFEMFunctionReturn(0);
1287}
1288
1289template <int SPACE_DIM, typename OpBase>
1290MoFEMErrorCode OpMixDivTimesScalarImpl<SPACE_DIM, GAUSS, OpBase>::iNtegrate(
1291 EntitiesFieldData::EntData &row_data,
1292 EntitiesFieldData::EntData &col_data) {
1293 MoFEMFunctionBegin;
1294
1295 auto t_w = this->getFTensor0IntegrationWeight();
1296
1297 size_t nb_base_functions = row_data.getN().size2() / 3;
1298 auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1299
1300 const double alpha_constant = alphaConstant();
1301
1302 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1303 const double alpha = alpha_constant * this->getMeasure() * t_w;
1304
1305 size_t rr = 0;
1306 for (; rr != OpBase::nbRows; ++rr) {
1307 const double t_row_div_base = t_row_diff_base(i, i);
1308 auto t_col_base = col_data.getFTensor0N(gg, 0);
1309 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1310 this->locMat(rr, cc) += alpha * t_row_div_base * t_col_base;
1311 ++t_col_base;
1312 }
1313 ++t_row_diff_base;
1314 }
1315 for (; rr < nb_base_functions; ++rr)
1316 ++t_row_diff_base;
1317
1318 ++t_w;
1319 }
1320
1321 MoFEMFunctionReturn(0);
1322}
1323
1324template <int SPACE_DIM, typename OpBase, CoordinateTypes CoordSys>
1325MoFEMErrorCode
1326OpMixDivTimesVecImpl<SPACE_DIM, GAUSS, OpBase, CoordSys>::iNtegrate(
1327 EntitiesFieldData::EntData &row_data,
1328 EntitiesFieldData::EntData &col_data) {
1329 MoFEMFunctionBegin;
1330
1331 auto t_w = this->getFTensor0IntegrationWeight();
1332 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1333
1334 size_t nb_base_functions = row_data.getN().size2() / 3;
1335 auto t_row_diff_base = row_data.getFTensor2DiffN<3, SPACE_DIM>();
1336 auto t_row_base = row_data.getFTensor1N<3>();
1337 const double alpha_constant = alphaConstant() * this->getMeasure();
1338 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1339
1340 const double alpha =
1341 alpha_constant * t_w * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1342
1343 size_t rr = 0;
1344 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1345 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1346 this->locMat, SPACE_DIM * rr);
1347 const double t_row_div_base = t_row_diff_base(i, i);
1348 auto t_col_base = col_data.getFTensor0N(gg, 0);
1349
1350 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1351 t_mat_diag(i) += alpha * t_row_div_base * t_col_base;
1352 if constexpr (CoordSys == CYLINDRICAL) {
1353 t_mat_diag(i) += t_row_base(0) * (alpha / t_coords(0)) * t_col_base;
1354 }
1355 ++t_col_base;
1356 ++t_mat_diag;
1357 }
1358 ++t_row_base;
1359 ++t_row_diff_base;
1360 }
1361 for (; rr < nb_base_functions; ++rr) {
1362 ++t_row_diff_base;
1363 ++t_row_base;
1364 }
1365
1366 ++t_w;
1367 ++t_coords;
1368 }
1369
1370 MoFEMFunctionReturn(0);
1371}
1372
1373template <int SPACE_DIM, typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
1374MoFEMErrorCode
1375OpMixScalarTimesDivImpl<SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM>::iNtegrate(
1376 EntitiesFieldData::EntData &row_data,
1377 EntitiesFieldData::EntData &col_data) {
1378 MoFEMFunctionBegin;
1379
1380#ifndef NDEBUG
1381 if (OpBase::locMat.size2() % SPACE_DIM)
1382 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1383 "Number of rows in matrix should be multiple of space dimensions");
1384#endif
1385
1386 // When we move to C++17 add if constexpr()
1387 if constexpr (COORDINATE_SYSTEM == POLAR || COORDINATE_SYSTEM == SPHERICAL)
1388 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1389 "%s coordiante not implemented",
1390 CoordinateTypesNames[COORDINATE_SYSTEM]);
1391
1392 auto t_w = this->getFTensor0IntegrationWeight();
1393 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1394 size_t nb_base_functions_row = row_data.getN().size2();
1395 auto t_row_base = row_data.getFTensor0N();
1396 const double vol = this->getMeasure();
1397 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1398
1399 const double alpha =
1400 alphaConstant(t_coords(0), t_coords(1), t_coords(2)) * t_w * vol;
1401
1402 size_t rr = 0;
1403 auto t_m = getFTensor1FromPtr<SPACE_DIM>(OpBase::locMat.data().data());
1404
1405 // When we move to C++17 add if constexpr()
1406 if constexpr (COORDINATE_SYSTEM == CARTESIAN) {
1407 for (; rr != OpBase::nbRows; ++rr) {
1408 const double r_val = alpha * t_row_base;
1409 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1410 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1411 t_m(i) += r_val * t_col_diff_base(i);
1412 ++t_col_diff_base;
1413 ++t_m;
1414 }
1415 ++t_row_base;
1416 }
1417 }
1418
1419 // When we move to C++17 add if constexpr()
1420 if constexpr (COORDINATE_SYSTEM == CYLINDRICAL) {
1421 for (; rr != OpBase::nbRows; ++rr) {
1422 const double r_val = alpha * t_row_base;
1423 auto t_col_base = col_data.getFTensor0N(gg, 0);
1424 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1425 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1426 t_m(i) += r_val * t_col_diff_base(i);
1427 t_m(0) += (r_val / t_coords(0)) * t_col_base;
1428 ++t_col_base;
1429 ++t_col_diff_base;
1430 ++t_m;
1431 }
1432 ++t_row_base;
1433 }
1434 }
1435
1436 for (; rr < nb_base_functions_row; ++rr)
1437 ++t_row_base;
1438
1439 ++t_w;
1440 ++t_coords;
1441 }
1442
1443 MoFEMFunctionReturn(0);
1444}
1445
1446template <int SPACE_DIM, typename OpBase>
1447MoFEMErrorCode
1448OpMixVectorTimesGradImpl<3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase>::iNtegrate(
1449 EntitiesFieldData::EntData &row_data,
1450 EntitiesFieldData::EntData &col_data) {
1451 MoFEMFunctionBegin;
1452
1453 auto t_w = this->getFTensor0IntegrationWeight();
1454
1455 size_t nb_base_functions = row_data.getN().size2() / 3;
1456 auto t_row_base = row_data.getFTensor1N<3>();
1457 auto &mat = this->locMat;
1458 const double alpha_constant = alphaConstant();
1459 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1460
1461 const double alpha = alpha_constant * this->getMeasure() * t_w;
1462
1463 size_t rr = 0;
1464 for (; rr != OpBase::nbRows; ++rr) {
1465 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1466 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1467 mat(rr, cc) += alpha * t_row_base(i) * t_col_diff_base(i);
1468 ++t_col_diff_base;
1469 }
1470 ++t_row_base;
1471 }
1472 for (; rr < nb_base_functions; ++rr)
1473 ++t_row_base;
1474
1475 ++t_w;
1476 }
1477
1478 MoFEMFunctionReturn(0);
1479}
1480
1481template <int SPACE_DIM, typename OpBase>
1482MoFEMErrorCode
1483OpMixVectorTimesGradImpl<1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase>::iNtegrate(
1484 EntitiesFieldData::EntData &row_data,
1485 EntitiesFieldData::EntData &col_data) {
1486 MoFEMFunctionBegin;
1487
1488 auto t_w = this->getFTensor0IntegrationWeight();
1489
1490 size_t nb_base_functions = row_data.getN().size2();
1491 auto t_row_base = row_data.getFTensor0N();
1492
1493 auto get_t_vec = [&](const int rr) {
1494 std::array<double *, SPACE_DIM> ptrs;
1495 for (auto i = 0; i != SPACE_DIM; ++i)
1496 ptrs[i] = &OpBase::locMat(rr + i, 0);
1497 return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, SPACE_DIM>(ptrs);
1498 };
1499
1500 const double alpha_constant = alphaConstant();
1501 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1502
1503 const double alpha = alpha_constant * this->getMeasure() * t_w;
1504
1505 size_t rr = 0;
1506 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1507 auto t_vec = get_t_vec(SPACE_DIM * rr);
1508 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1509 for (size_t cc = 0; cc != OpBase::nbCols; ++cc) {
1510 t_vec(i) += alpha * t_row_base * t_col_diff_base(i);
1511 ++t_col_diff_base;
1512 ++t_vec;
1513 }
1514 ++t_row_base;
1515 }
1516 for (; rr < nb_base_functions; ++rr)
1517 ++t_row_base;
1518
1519 ++t_w;
1520 }
1521
1522 MoFEMFunctionReturn(0);
1523}
1524
1525template <int SPACE_DIM, typename OpBase>
1526MoFEMErrorCode OpMixTensorTimesGradImpl<SPACE_DIM, GAUSS, OpBase>::iNtegrate(
1527 EntitiesFieldData::EntData &row_data,
1528 EntitiesFieldData::EntData &col_data) {
1529 MoFEMFunctionBegin;
1530
1531 auto t_w = this->getFTensor0IntegrationWeight();
1532 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1533
1534 size_t nb_base_functions = row_data.getN().size2() / 3;
1535 auto t_row_base = row_data.getFTensor1N<3>();
1536 const double alpha_constant = alphaConstant() * this->getMeasure();
1537 for (size_t gg = 0; gg != OpBase::nbIntegrationPts; ++gg) {
1538
1539 const double alpha =
1540 alpha_constant * betaCoeff(t_coords(0), t_coords(1), t_coords(2)) * t_w;
1541
1542 size_t rr = 0;
1543 for (; rr != OpBase::nbRows / SPACE_DIM; ++rr) {
1544 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1545 this->locMat, SPACE_DIM * rr);
1546 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1547
1548 for (size_t cc = 0; cc != OpBase::nbCols / SPACE_DIM; ++cc) {
1549 t_mat_diag(i) += alpha * t_row_base(j) * t_col_diff_base(j);
1550 ++t_col_diff_base;
1551 ++t_mat_diag;
1552 }
1553
1554 ++t_row_base;
1555 }
1556 for (; rr < nb_base_functions; ++rr)
1557 ++t_row_base;
1558
1559 ++t_w;
1560 ++t_coords;
1561 }
1562
1563 MoFEMFunctionReturn(0);
1564}
1565
1566} // namespace MoFEM
1567
1568#endif //__BILINEAR_FORMS_INTEGRATORS_IMPL_HPP__
a
constexpr double a
Definition approx_sphere.cpp:30
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
FIELD_DIM
constexpr int FIELD_DIM
Definition child_and_parent.cpp:16
FTensor::Kronecker_Delta
Kronecker Delta class.
Definition Kronecker_Delta.hpp:15
FTensor::Tensor2_symmetric
Definition Tensor2_symmetric_value.hpp:14
CoordinateTypesNames
static const char *const CoordinateTypesNames[]
Coordinate system names.
Definition definitions.h:121
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
CoordinateTypes
CoordinateTypes
Coodinate system.
Definition definitions.h:127
CYLINDRICAL
@ CYLINDRICAL
Definition definitions.h:130
POLAR
@ POLAR
Definition definitions.h:129
CARTESIAN
@ CARTESIAN
Definition definitions.h:128
SPHERICAL
@ SPHERICAL
Definition definitions.h:131
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
BASE_DIM
constexpr int BASE_DIM
Definition dg_projection.cpp:35
t_kd
constexpr auto t_kd
Definition free_surface.cpp:161
MoFEM::IntegrationType
IntegrationType
Form integrator integration types.
Definition FormsIntegrators.hpp:237
MoFEM::ScalarFun
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
Definition FormsIntegrators.hpp:249
MoFEM::scalar_fun_one
double scalar_fun_one(const double, const double, const double)
Default scalar function returning 1.
Definition FormsIntegrators.hpp:260
MoFEM::GAUSS
@ GAUSS
Gaussian quadrature integration.
Definition FormsIntegrators.hpp:238
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
J
FTensor::Index< 'J', DIM1 > J
Definition level_set.cpp:30
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:23
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::GetFTensor4FromMatType
decltype(GetFTensor4FromMatImpl< Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3, S, DL, M >::get(std::declval< M & >(), 0, 0)) GetFTensor4FromMatType
Definition Templates.hpp:1360
MoFEM::CacheMatsTypeType
std::map< std::pair< int, int >, boost::shared_ptr< MatrixDouble > > CacheMatsTypeType
Definition BiLinearFormsIntegratorsImpl.hpp:19
MoFEM::ConstantFun
boost::function< double()> ConstantFun
Constant function type.
Definition FormsIntegrators.hpp:280
I
constexpr IntegrationType I
Definition operators_tests.cpp:31
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition EntitiesFieldData.hpp:150
MoFEM::EntitiesFieldData::EntData::getFTensor2DiffN
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
Definition EntitiesFieldData.cpp:297
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition EntitiesFieldData.hpp:1806
MoFEM::EntitiesFieldData::EntData::getFieldEntities
const VectorFieldEntities & getFieldEntities() const
Get field entities (const version)
Definition EntitiesFieldData.hpp:1591
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
auto getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition EntitiesFieldData.hpp:802
MoFEM::EntitiesFieldData::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition EntitiesFieldData.hpp:1600
MoFEM::EntitiesFieldData::EntData::getFTensor1N
auto getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
Definition EntitiesFieldData.hpp:948
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Definition EntitiesFieldData.hpp:1523
MoFEM::OpBaseImpl::rowSide
int rowSide
row side number
Definition FormsIntegrators.hpp:373
MoFEM::OpBaseImpl::colSide
int colSide
column side number
Definition FormsIntegrators.hpp:374
MoFEM::OpBaseImpl::nbRows
int nbRows
number of dofs on rows
Definition FormsIntegrators.hpp:368
MoFEM::OpBaseImpl::colType
EntityType colType
column type
Definition FormsIntegrators.hpp:376
MoFEM::OpBaseImpl::nbIntegrationPts
int nbIntegrationPts
number of integration points
Definition FormsIntegrators.hpp:370
MoFEM::OpBaseImpl::locMat
MatrixDouble locMat
local entity block matrix
Definition FormsIntegrators.hpp:381
MoFEM::OpBaseImpl::nbCols
int nbCols
number if dof on column
Definition FormsIntegrators.hpp:369
MoFEM::OpBaseImpl::nbRowBaseFunctions
int nbRowBaseFunctions
number or row base functions
Definition FormsIntegrators.hpp:371
MoFEM::OpBaseImpl::rowType
EntityType rowType
row type
Definition FormsIntegrators.hpp:375
MoFEM::OpGradGradImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::OpGradGradImpl
OpGradGradImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition BiLinearFormsIntegratorsImpl.hpp:28
MoFEM::OpGradGradImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:27
MoFEM::OpGradGradImpl< 1, 1, SPACE_DIM, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition BiLinearFormsIntegratorsImpl.hpp:39
MoFEM::OpGradGradImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition BiLinearFormsIntegratorsImpl.hpp:58
MoFEM::OpGradGradImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::OpGradGradImpl
OpGradGradImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition BiLinearFormsIntegratorsImpl.hpp:47
MoFEM::OpGradGradImpl< 1, FIELD_DIM, SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:46
MoFEM::OpGradGradSymTensorGradGradImpl< 1, FIELD_DIM, SPACE_DIM, S, GAUSS, OpBase >::matD
boost::shared_ptr< MatrixDouble > matD
Definition BiLinearFormsIntegratorsImpl.hpp:252
MoFEM::OpGradGradSymTensorGradGradImpl< 1, FIELD_DIM, SPACE_DIM, S, GAUSS, OpBase >::OpGradGradSymTensorGradGradImpl
OpGradGradSymTensorGradGradImpl(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mat_D, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition BiLinearFormsIntegratorsImpl.hpp:241
MoFEM::OpGradGradTensorGradGradImpl< 1, FIELD_DIM, SPACE_DIM, S, GAUSS, OpBase >::matD
boost::shared_ptr< MatrixDouble > matD
Definition BiLinearFormsIntegratorsImpl.hpp:275
MoFEM::OpGradGradTensorGradGradImpl< 1, FIELD_DIM, SPACE_DIM, S, GAUSS, OpBase >::OpGradGradTensorGradGradImpl
OpGradGradTensorGradGradImpl(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mat_D, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition BiLinearFormsIntegratorsImpl.hpp:264
MoFEM::OpGradSymTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition BiLinearFormsIntegratorsImpl.hpp:207
MoFEM::OpGradSymTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:193
MoFEM::OpGradSymTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::matD
boost::shared_ptr< MatrixDouble > matD
Definition BiLinearFormsIntegratorsImpl.hpp:206
MoFEM::OpGradSymTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::OpGradSymTensorGradImpl
OpGradSymTensorGradImpl(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mat_D, boost::shared_ptr< Range > ents_ptr=nullptr, ScalarFun beta=scalar_fun_one)
Definition BiLinearFormsIntegratorsImpl.hpp:194
MoFEM::OpGradTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::matD
boost::shared_ptr< MatrixDouble > matD
Definition BiLinearFormsIntegratorsImpl.hpp:228
MoFEM::OpGradTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition BiLinearFormsIntegratorsImpl.hpp:229
MoFEM::OpGradTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:219
MoFEM::OpGradTensorGradImpl< 1, SPACE_DIM, SPACE_DIM, S, GAUSS, OpBase >::OpGradTensorGradImpl
OpGradTensorGradImpl(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< MatrixDouble > mat_D, ScalarFun beta=scalar_fun_one)
Definition BiLinearFormsIntegratorsImpl.hpp:220
MoFEM::OpMassCacheImpl< 1, FIELD_DIM, I, OpBase >::scalarBeta
double scalarBeta
Definition BiLinearFormsIntegratorsImpl.hpp:182
MoFEM::OpMassCacheImpl< 1, FIELD_DIM, I, OpBase >::cacheLocMats
CacheMatsTypeType cacheLocMats
Definition BiLinearFormsIntegratorsImpl.hpp:183
MoFEM::OpMassCacheImpl< 1, FIELD_DIM, I, OpBase >::OpMassCacheImpl
OpMassCacheImpl(CacheMatsTypeType cache, const std::string row_field_name, const std::string col_field_name, const double beta, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition BiLinearFormsIntegratorsImpl.hpp:170
MoFEM::OpMassImpl< 1, 1, GAUSS, OpBase >::OpMassImpl
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=[](double, double, double) constexpr { return 1;}, boost::shared_ptr< Range > ents_ptr=nullptr, boost::shared_ptr< MatrixDouble > cache_mat=nullptr)
Definition BiLinearFormsIntegratorsImpl.hpp:69
MoFEM::OpMassImpl< 1, 1, GAUSS, OpBase >::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
Definition BiLinearFormsIntegratorsImpl.hpp:87
MoFEM::OpMassImpl< 1, FIELD_DIM, GAUSS, OpBase >::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
Class dedicated to integrate operator.
Definition BiLinearFormsIntegratorsImpl.hpp:107
MoFEM::OpMassImpl< 1, FIELD_DIM, GAUSS, OpBase >::iNtegrate
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
Definition BiLinearFormsIntegratorsImpl.hpp:103
MoFEM::OpMassImpl< 3, 4, GAUSS, OpBase >::OpMassImpl
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition BiLinearFormsIntegratorsImpl.hpp:132
MoFEM::OpMassImpl< 3, 9, GAUSS, OpBase >::OpMassImpl
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition BiLinearFormsIntegratorsImpl.hpp:149
MoFEM::OpMassImpl< 3, FIELD_DIM, GAUSS, OpBase >::betaCoeff
ScalarFun betaCoeff
Definition BiLinearFormsIntegratorsImpl.hpp:125
MoFEM::OpMassImpl< 3, FIELD_DIM, GAUSS, OpBase >::OpMassImpl
OpMassImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun beta=scalar_fun_one, boost::shared_ptr< Range > ents_ptr=nullptr)
Definition BiLinearFormsIntegratorsImpl.hpp:115
MoFEM::OpMixDivTimesScalarImpl< SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition BiLinearFormsIntegratorsImpl.hpp:298
MoFEM::OpMixDivTimesScalarImpl< SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:297
MoFEM::OpMixDivTimesScalarImpl< SPACE_DIM, GAUSS, OpBase >::OpMixDivTimesScalarImpl
OpMixDivTimesScalarImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose=false, const bool only_transpose=false)
Definition BiLinearFormsIntegratorsImpl.hpp:285
MoFEM::OpMixDivTimesVecImpl< SPACE_DIM, GAUSS, OpBase, CoordSys >::OpMixDivTimesVecImpl
OpMixDivTimesVecImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose, const bool only_transpose=false)
Definition BiLinearFormsIntegratorsImpl.hpp:311
MoFEM::OpMixDivTimesVecImpl< SPACE_DIM, GAUSS, OpBase, CoordSys >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:332
MoFEM::OpMixDivTimesVecImpl< SPACE_DIM, GAUSS, OpBase, CoordSys >::OpMixDivTimesVecImpl
OpMixDivTimesVecImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, ScalarFun beta_fun, const bool assemble_transpose, const bool only_transpose=false)
Definition BiLinearFormsIntegratorsImpl.hpp:321
MoFEM::OpMixScalarTimesDivImpl< SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM >::alphaConstant
ScalarFun alphaConstant
Definition BiLinearFormsIntegratorsImpl.hpp:361
MoFEM::OpMixScalarTimesDivImpl< SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM >::OpMixScalarTimesDivImpl
OpMixScalarTimesDivImpl(const std::string row_field_name, const std::string col_field_name, ScalarFun alpha_fun=[](double, double, double) constexpr { return 1;}, const bool assemble_transpose=false, const bool only_transpose=false)
Definition BiLinearFormsIntegratorsImpl.hpp:348
MoFEM::OpMixScalarTimesDivImpl< SPACE_DIM, GAUSS, OpBase, COORDINATE_SYSTEM >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:360
MoFEM::OpMixTensorTimesGradImpl< SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:440
MoFEM::OpMixTensorTimesGradImpl< SPACE_DIM, GAUSS, OpBase >::OpMixTensorTimesGradImpl
OpMixTensorTimesGradImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose, const bool only_transpose=false)
Definition BiLinearFormsIntegratorsImpl.hpp:418
MoFEM::OpMixTensorTimesGradImpl< SPACE_DIM, GAUSS, OpBase >::OpMixTensorTimesGradImpl
OpMixTensorTimesGradImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, ScalarFun beta_coeff, const bool assemble_transpose, const bool only_transpose=false)
Definition BiLinearFormsIntegratorsImpl.hpp:428
MoFEM::OpMixTensorTimesGradImpl< SPACE_DIM, GAUSS, OpBase >::j
FTensor::Index< 'j', SPACE_DIM > j
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:441
MoFEM::OpMixVectorTimesGradImpl< 1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition BiLinearFormsIntegratorsImpl.hpp:407
MoFEM::OpMixVectorTimesGradImpl< 1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:406
MoFEM::OpMixVectorTimesGradImpl< 1, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::OpMixVectorTimesGradImpl
OpMixVectorTimesGradImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose=false, const bool only_transpose=false)
Definition BiLinearFormsIntegratorsImpl.hpp:394
MoFEM::OpMixVectorTimesGradImpl< 3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::i
FTensor::Index< 'i', SPACE_DIM > i
summit Index
Definition BiLinearFormsIntegratorsImpl.hpp:385
MoFEM::OpMixVectorTimesGradImpl< 3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::alphaConstant
ConstantFun alphaConstant
Definition BiLinearFormsIntegratorsImpl.hpp:386
MoFEM::OpMixVectorTimesGradImpl< 3, SPACE_DIM, SPACE_DIM, GAUSS, OpBase >::OpMixVectorTimesGradImpl
OpMixVectorTimesGradImpl(const std::string row_field_name, const std::string col_field_name, ConstantFun alpha_fun, const bool assemble_transpose=false, const bool only_transpose=false)
Definition BiLinearFormsIntegratorsImpl.hpp:373
Generated by Doxygen 1.9.8 and hosted at University of Glasgow