13 #ifndef __BILINEAR_FORMS_INTEGRATORS_HPP__
14 #define __BILINEAR_FORMS_INTEGRATORS_HPP__
22 template <
int SPACE_DIM,
typename OpBase>
26 const std::string col_field_name,
28 boost::shared_ptr<Range> ents_ptr =
nullptr)
29 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL, ents_ptr),
31 if (row_field_name == col_field_name)
41 template <
int FIELD_DIM,
int SPACE_DIM,
typename OpBase>
45 const std::string col_field_name,
47 boost::shared_ptr<Range> ents_ptr =
nullptr)
48 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL, ents_ptr),
50 if (row_field_name == col_field_name)
60 template <
int BASE_DIM,
int FIELD_DIM, IntegrationType I,
typename OpBase>
63 template <
typename OpBase>
67 const std::string row_field_name,
const std::string col_field_name,
68 ScalarFun beta = [](
double,
double,
double) constexpr {
return 1; },
69 boost::shared_ptr<Range> ents_ptr =
nullptr)
70 :
OpBase(row_field_name, col_field_name, OpBase::OPROWCOL, ents_ptr),
72 if (row_field_name == col_field_name)
82 template <
int FIELD_DIM,
typename OpBase>
93 template <
int FIELD_DIM,
typename OpBase>
96 OpMassImpl(
const std::string row_field_name,
const std::string col_field_name,
98 boost::shared_ptr<Range> ents_ptr =
nullptr)
99 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL, ents_ptr),
101 if (row_field_name == col_field_name)
111 template <
typename OpBase>
113 OpMassImpl(
const std::string row_field_name,
const std::string col_field_name,
115 boost::shared_ptr<Range> ents_ptr =
nullptr)
116 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL, ents_ptr),
118 if (row_field_name == col_field_name)
128 template <
typename OpBase>
130 OpMassImpl(
const std::string row_field_name,
const std::string col_field_name,
132 boost::shared_ptr<Range> ents_ptr =
nullptr)
133 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL, ents_ptr),
135 if (row_field_name == col_field_name)
149 template <
int SPACE_DIM,
int S,
typename OpBase>
154 const std::string col_field_name,
155 boost::shared_ptr<MatrixDouble> mat_D,
156 boost::shared_ptr<Range> ents_ptr =
nullptr,
158 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL, ents_ptr),
159 matD(mat_D), betaCoeff(beta) {
160 if (row_field_name == col_field_name)
165 boost::shared_ptr<MatrixDouble>
matD;
175 template <
int SPACE_DIM,
int S,
typename OpBase>
180 const std::string col_field_name,
181 boost::shared_ptr<MatrixDouble> mat_D,
183 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL), matD(mat_D),
187 boost::shared_ptr<MatrixDouble>
matD;
197 template <
int SPACE_DIM,
int S,
typename OpBase>
202 const std::string col_field_name,
203 boost::shared_ptr<MatrixDouble> mat_D,
204 boost::shared_ptr<Range> ents_ptr =
nullptr)
205 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL, ents_ptr),
207 if (row_field_name == col_field_name)
212 boost::shared_ptr<MatrixDouble>
matD;
217 template <
int SPACE_DIM, IntegrationType I,
typename OpBase>
220 template <
int SPACE_DIM,
typename OpBase>
223 const std::string col_field_name,
225 const bool assemble_transpose =
false,
226 const bool only_transpose =
false)
227 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL),
228 alphaConstant(alpha_fun) {
229 this->assembleTranspose = assemble_transpose;
230 this->onlyTranspose = only_transpose;
244 template <
int SPACE_DIM,
typename OpBase, CoordinateTypes CoordSys>
249 const std::string col_field_name,
ConstantFun alpha_fun,
250 const bool assemble_transpose,
251 const bool only_transpose =
false)
252 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL),
253 alphaConstant(alpha_fun) {
254 this->assembleTranspose = assemble_transpose;
255 this->onlyTranspose = only_transpose;
259 const std::string col_field_name,
ConstantFun alpha_fun,
260 ScalarFun beta_fun,
const bool assemble_transpose,
261 const bool only_transpose =
false)
262 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL),
263 alphaConstant(alpha_fun), betaCoeff(beta_fun) {
264 this->assembleTranspose = assemble_transpose;
265 this->onlyTranspose = only_transpose;
282 template <
int SPACE_DIM,
typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
286 const std::string row_field_name,
const std::string col_field_name,
287 ScalarFun alpha_fun = [](
double,
double,
double) constexpr {
return 1; },
288 const bool assemble_transpose =
false,
const bool only_transpose =
false)
289 :
OpBase(row_field_name, col_field_name, OpBase::OPROWCOL),
290 alphaConstant(alpha_fun) {
291 this->assembleTranspose = assemble_transpose;
292 this->onlyTranspose = only_transpose;
307 template <
int SPACE_DIM,
typename OpBase>
311 const std::string col_field_name,
313 const bool assemble_transpose =
false,
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;
328 template <
int SPACE_DIM,
typename OpBase>
332 const std::string col_field_name,
334 const bool assemble_transpose =
false,
335 const bool only_transpose =
false)
336 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL),
337 alphaConstant(alpha_fun) {
338 this->assembleTranspose = assemble_transpose;
339 this->onlyTranspose = only_transpose;
349 template <
int SPACE_DIM, IntegrationType I,
typename OpBase>
352 template <
int SPACE_DIM,
typename OpBase>
356 const std::string col_field_name,
357 ConstantFun alpha_fun,
const bool assemble_transpose,
358 const bool only_transpose =
false)
359 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL),
360 alphaConstant(alpha_fun) {
361 this->assembleTranspose = assemble_transpose;
362 this->onlyTranspose = only_transpose;
366 const std::string col_field_name,
368 const bool assemble_transpose,
369 const bool only_transpose =
false)
370 :
OpBase(row_field_name, col_field_name,
OpBase::OPROWCOL),
371 alphaConstant(alpha_fun), betaCoeff(beta_coeff) {
372 this->assembleTranspose = assemble_transpose;
373 this->onlyTranspose = only_transpose;
393 template <
typename EleOp>
394 template <AssemblyType A>
395 template <IntegrationType I>
404 template <
int BASE_DIM,
int FIELD_DIM,
int SPACE_DIM>
420 template <
int BASE_DIM,
int FIELD_DIM>
433 template <
int BASE_DIM,
int FIELD_DIM,
int SPACE_DIM,
int S = 0>
447 template <
int BASE_DIM,
int FIELD_DIM,
int SPACE_DIM,
int S = 0>
462 template <
int BASE_DIM,
int FIELD_DIM,
int SPACE_DIM,
int S = 0>
471 template <
int SPACE_DIM>
479 template <
int SPACE_DIM, CoordinateTypes CoordSys = CARTESIAN>
487 template <
int SPACE_DIM, CoordinateTypes COORDINATE_SYSTEM = CARTESIAN>
496 template <
int BASE_DIM,
int FIELD_DIM,
int SPACE_DIM>
505 template <
int SPACE_DIM>
509 template <
int SPACE_DIM,
typename OpBase>
516 const double vol = OpBase::getMeasure();
518 auto t_w = OpBase::getFTensor0IntegrationWeight();
522 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
526 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
528 const double alpha = t_w * beta;
553 template <
int FIELD_DIM,
int SPACE_DIM,
typename OpBase>
561 const double vol = OpBase::getMeasure();
563 auto t_w = OpBase::getFTensor0IntegrationWeight();
567 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
572 auto get_t_vec = [&](
const int rr) {
573 std::array<double *, FIELD_DIM> ptrs;
582 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
584 const double alpha = t_w * beta;
595 t_vec(
j) += alpha * (t_row_grad(
i) * t_col_grad(
i));
612 template <
typename OpBase>
619 auto log_error = [&]() {
629 "Wrong number of base functions on rows %d < %d",
635 "Wrong number of base functions on cols %d < %d",
641 "Wrong number of integration points on rows %d != %d",
647 "Wrong number of integration points on cols %d != %d",
653 const double vol = OpBase::getMeasure();
655 auto t_w = OpBase::getFTensor0IntegrationWeight();
659 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
662 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
664 const double alpha = t_w * beta;
674 *a_mat_ptr += alpha * (t_row_base * t_col_base);
688 template <
int FIELD_DIM,
typename OpBase>
694 const double vol = OpBase::getMeasure();
696 auto t_w = OpBase::getFTensor0IntegrationWeight();
700 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
703 auto get_t_vec = [&](
const int rr) {
704 std::array<double *, FIELD_DIM> ptrs;
713 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
715 const double alpha = t_w * beta;
726 t_vec(
i) += alpha * (t_row_base * t_col_base);
740 template <
int FIELD_DIM,
typename OpBase>
746 size_t nb_base_functions = row_data.
getN().size2() / 3;
748 const double vol = OpBase::getMeasure();
750 auto t_w = OpBase::getFTensor0IntegrationWeight();
754 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
757 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
759 const double alpha = t_w * beta;
769 (*a_mat_ptr) += alpha * (t_row_base(
i) * t_col_base(
i));
775 for (; rr < nb_base_functions; ++rr)
783 template <
typename OpBase>
790 auto get_t_vec = [&](
const int rr) {
794 size_t nb_base_functions = row_data.
getN().size2() / 3;
796 const double vol = OpBase::getMeasure();
798 auto t_w = OpBase::getFTensor0IntegrationWeight();
802 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
805 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
807 const double alpha = t_w * beta;
813 auto t_vec = get_t_vec(2 * rr);
817 t_vec(
I) += alpha * (t_row_base(
k) * t_col_base(
k));
823 for (; rr < nb_base_functions; ++rr)
831 template <
typename OpBase>
838 auto get_t_vec = [&](
const int rr) {
843 size_t nb_base_functions = row_data.
getN().size2() / 3;
845 const double vol = OpBase::getMeasure();
847 auto t_w = OpBase::getFTensor0IntegrationWeight();
851 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
854 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
856 const double alpha = t_w * beta;
862 auto t_vec = get_t_vec(3 * rr);
866 t_vec(
i) += alpha * (t_row_base(
k) * t_col_base(
k));
872 for (; rr < nb_base_functions; ++rr)
880 template <
int SPACE_DIM,
int S,
typename OpBase>
887 const size_t nb_row_dofs = row_data.
getIndices().size();
888 const size_t nb_col_dofs = col_data.
getIndices().size();
890 if (nb_row_dofs && nb_col_dofs) {
901 double vol = OpBase::getMeasure();
904 auto t_w = OpBase::getFTensor0IntegrationWeight();
911 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
914 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
920 double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
927 auto t_m = OpBase::template getLocMat<SPACE_DIM>(
SPACE_DIM * rr);
933 t_rowD(
l,
j,
k) = t_D(
i,
j,
k,
l) * (
a * t_row_diff_base(
i));
940 for (
int cc = 0; cc <= nb_cols; ++cc) {
943 t_m(
i,
j) += t_rowD(
i,
j,
k) * t_col_diff_base(
k);
968 auto t_m_rr = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
970 auto t_m_cc = getFTensor2FromArray<SPACE_DIM, SPACE_DIM>(
974 t_m_rr(
i,
j) = t_m_cc(
j,
i);
985 template <
int SPACE_DIM,
int S,
typename OpBase>
999 auto &row_hessian = row_data.
getN(BaseDerivatives::SecondDerivative);
1000 auto &col_hessian = col_data.
getN(BaseDerivatives::SecondDerivative);
1005 "Wrong number of integration pts (%d != %d)",
1008 if (row_hessian.size2() !=
1011 "Wrong number of base functions (%d != %d)",
1017 "Wrong number of base functions (%d < %d)", row_hessian.size2(),
1022 "Wrong number of integration pts (%d != %d)",
1027 "Wrong number of base functions (%d < %d)", col_hessian.size2(),
1033 double vol = OpBase::getMeasure();
1036 auto t_w = OpBase::getFTensor0IntegrationWeight();
1038 auto t_row_diff2 = getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(
1039 &*row_hessian.data().begin());
1043 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, S>(*matD);
1049 double a = t_w * vol;
1059 t_rowD(
k,
l) = t_D(
i,
j,
k,
l) * (
a * t_row_diff2(
i,
j));
1063 getFTensor2SymmetricLowerFromPtr<SPACE_DIM>(&col_hessian(gg, 0));
1069 *m_ptr += t_rowD(
i,
j) * t_col_diff2(
i,
j);
1094 template <
int SPACE_DIM,
int S,
typename OpBase>
1101 const size_t nb_row_dofs = row_data.
getIndices().size();
1102 const size_t nb_col_dofs = col_data.
getIndices().size();
1104 if (nb_row_dofs && nb_col_dofs) {
1112 double vol = OpBase::getMeasure();
1115 auto t_w = OpBase::getFTensor0IntegrationWeight();
1122 getFTensor4FromMat<SPACE_DIM, SPACE_DIM, SPACE_DIM, SPACE_DIM, S>(
1126 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1132 double a = t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1139 auto t_m = OpBase::template getLocMat<SPACE_DIM>(
SPACE_DIM * rr);
1143 t_row(
i,
k,
l) = t_D(
i,
j,
k,
l) * (
a * t_row_diff_base(
j));
1152 t_m(
i,
k) += t_row(
i,
k,
l) * t_col_diff_base(
l);
1178 template <
int SPACE_DIM,
typename OpBase>
1184 auto t_w = this->getFTensor0IntegrationWeight();
1186 size_t nb_base_functions = row_data.
getN().size2() / 3;
1189 const double alpha_constant = alphaConstant();
1192 const double alpha = alpha_constant * this->getMeasure() * t_w;
1196 const double t_row_div_base = t_row_diff_base(
i,
i);
1199 this->locMat(rr, cc) += alpha * t_row_div_base * t_col_base;
1204 for (; rr < nb_base_functions; ++rr)
1213 template <
int SPACE_DIM,
typename OpBase, CoordinateTypes CoordSys>
1220 auto t_w = this->getFTensor0IntegrationWeight();
1221 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1223 size_t nb_base_functions = row_data.
getN().size2() / 3;
1226 const double alpha_constant = alphaConstant() * this->getMeasure();
1229 const double alpha =
1230 alpha_constant * t_w * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1234 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1236 const double t_row_div_base = t_row_diff_base(
i,
i);
1240 t_mat_diag(
i) += alpha * t_row_div_base * t_col_base;
1242 t_mat_diag(
i) += t_row_base(0) * (alpha / t_coords(0)) * t_col_base;
1250 for (; rr < nb_base_functions; ++rr) {
1262 template <
int SPACE_DIM,
typename OpBase, CoordinateTypes COORDINATE_SYSTEM>
1272 "Number of rows in matrix should be multiple of space dimensions");
1276 if constexpr (COORDINATE_SYSTEM ==
POLAR || COORDINATE_SYSTEM ==
SPHERICAL)
1278 "%s coordiante not implemented",
1281 auto t_w = this->getFTensor0IntegrationWeight();
1282 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1283 size_t nb_base_functions_row = row_data.
getN().size2();
1285 const double vol = this->getMeasure();
1288 const double alpha =
1289 alphaConstant(t_coords(0), t_coords(1), t_coords(2)) * t_w * vol;
1292 auto t_m = getFTensor1FromPtr<SPACE_DIM>(
OpBase::locMat.data().data());
1295 if constexpr (COORDINATE_SYSTEM ==
CARTESIAN) {
1297 const double r_val = alpha * t_row_base;
1300 t_m(
i) += r_val * t_col_diff_base(
i);
1311 const double r_val = alpha * t_row_base;
1315 t_m(
i) += r_val * t_col_diff_base(
i);
1316 t_m(0) += (r_val / t_coords(0)) * t_col_base;
1325 for (; rr < nb_base_functions_row; ++rr)
1335 template <
int SPACE_DIM,
typename OpBase>
1342 auto t_w = this->getFTensor0IntegrationWeight();
1344 size_t nb_base_functions = row_data.
getN().size2() / 3;
1346 auto &mat = this->locMat;
1347 const double alpha_constant = alphaConstant();
1350 const double alpha = alpha_constant * this->getMeasure() * t_w;
1356 mat(rr, cc) += alpha * t_row_base(
i) * t_col_diff_base(
i);
1361 for (; rr < nb_base_functions; ++rr)
1370 template <
int SPACE_DIM,
typename OpBase>
1377 auto t_w = this->getFTensor0IntegrationWeight();
1379 size_t nb_base_functions = row_data.
getN().size2();
1382 auto get_t_vec = [&](
const int rr) {
1383 std::array<double *, SPACE_DIM> ptrs;
1389 const double alpha_constant = alphaConstant();
1392 const double alpha = alpha_constant * this->getMeasure() * t_w;
1399 t_vec(
i) += alpha * t_row_base * t_col_diff_base(
i);
1405 for (; rr < nb_base_functions; ++rr)
1414 template <
int SPACE_DIM,
typename OpBase>
1420 auto t_w = this->getFTensor0IntegrationWeight();
1421 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1423 size_t nb_base_functions = row_data.
getN().size2() / 3;
1425 const double alpha_constant = alphaConstant() * this->getMeasure();
1428 const double alpha =
1429 alpha_constant * betaCoeff(t_coords(0), t_coords(1), t_coords(2)) * t_w;
1433 auto t_mat_diag = getFTensor1FromArrayDiag<SPACE_DIM, SPACE_DIM>(
1438 t_mat_diag(
i) += alpha * t_row_base(
j) * t_col_diff_base(
j);
1445 for (; rr < nb_base_functions; ++rr)
1457 #endif //__BILINEAR_FORMS_INTEGRATORS_HPP__