7 #ifndef __LINEAR_FORMS_INTEGRATORS_HPP__
8 #define __LINEAR_FORMS_INTEGRATORS_HPP__
13 template <
typename OpBase>
struct S {
S() =
delete; };
18 template <
typename OpBase>
struct S {
S() =
delete; };
22 template <
int BASE_DIM,
int FIELD_DIM, IntegrationType I,
typename OpBase>
30 template <
typename OpBase>
42 boost::shared_ptr<Range> ents_ptr =
nullptr)
44 sourceFun(source_fun) {}
56 boost::shared_ptr<Range> ents_ptr =
nullptr)
58 sourceFun(source_fun) {}
65 template <
int FIELD_DIM,
typename OpBase>
79 boost::shared_ptr<Range> ents_ptr =
nullptr)
90 boost::shared_ptr<Range> ents_ptr =
nullptr)
92 sourceFun(source_fun) {}
99 template <
int FIELD_DIM,
typename OpBase>
105 boost::shared_ptr<Range> ents_ptr =
nullptr)
107 sourceFun(source_fun) {}
110 boost::shared_ptr<Range> ents_ptr =
nullptr)
112 sourceFun(source_fun) {}
119 template <
int BASE_DIM,
int S, IntegrationType I,
typename OpBase>
122 template <
int S,
typename OpBase>
126 const std::string
field_name, boost::shared_ptr<VectorDouble> vec,
127 ScalarFun beta_coeff = [](
double,
double,
double) constexpr {
return 1; },
128 boost::shared_ptr<Range> ents_ptr =
nullptr)
130 betaCoeff(beta_coeff) {}
142 template <
int FIELD_DIM,
int S,
typename OpBase>
146 const std::string
field_name, boost::shared_ptr<MatrixDouble> vec,
147 ScalarFun beta_coeff = [](
double,
double,
double) constexpr {
return 1; },
148 boost::shared_ptr<Range> ents_ptr =
nullptr)
150 betaCoeff(beta_coeff) {}
159 template <
int FIELD_DIM,
int S,
typename OpBase>
163 const std::string
field_name, boost::shared_ptr<MatrixDouble> vec,
164 ScalarFun beta_coeff = [](
double,
double,
double) constexpr {
return 1; },
165 boost::shared_ptr<Range> ents_ptr =
nullptr)
167 betaCoeff(beta_coeff) {}
180 template <
int SPACE_DIM,
int S,
typename OpBase>
187 const std::string
field_name, boost::shared_ptr<MatrixDouble> mat_vals,
188 ScalarFun beta_coeff = [](
double,
double,
double) constexpr {
return 1; },
189 boost::shared_ptr<Range> ents_ptr =
nullptr)
191 matVals(mat_vals), betaCoeff(beta_coeff) {}
199 template <
int SPACE_DIM,
int S,
typename OpBase>
207 const std::string
field_name, boost::shared_ptr<MatrixDouble> mat_vals,
208 ScalarFun beta_coeff = [](
double,
double,
double) constexpr {
return 1; },
209 boost::shared_ptr<Range> ents_ptr =
nullptr)
211 matVals(mat_vals), betaCoeff(beta_coeff) {}
223 template <
int SPACE_DIM,
int S,
typename OpBase>
228 const std::string
field_name, boost::shared_ptr<MatrixDouble> mat_vals,
229 ScalarFun beta_coeff = [](
double,
double,
double) constexpr {
return 1; })
231 betaCoeff(beta_coeff) {}
250 const std::string
field_name, boost::shared_ptr<MatrixDouble> mat_vals,
251 ScalarFun beta = [](
double,
double,
double) {
return 1; },
252 boost::shared_ptr<Range> ents_ptr =
nullptr)
254 matVals(mat_vals), betaConst(beta) {}
264 template <
int SPACE_DIM,
typename OpBase, CoordinateTypes CoordSys>
268 const std::string
field_name, boost::shared_ptr<VectorDouble> vec_vals,
269 ScalarFun beta = [](
double,
double,
double) constexpr {
return 1; },
270 boost::shared_ptr<Range> ents_ptr =
nullptr)
272 vecVals(vec_vals), betaConst(beta) {}
281 template <
int FIELD_DIM,
typename OpBase, CoordinateTypes CoordSys>
286 const std::string
field_name, boost::shared_ptr<VectorDouble> vec,
287 ScalarFun beta = [](
double,
double,
double) constexpr {
return 1; },
288 boost::shared_ptr<Range> ents_ptr =
nullptr)
309 template <
int SPACE_DIM, IntegrationType I,
typename OpBase>
312 template <
int SPACE_DIM,
typename OpBase>
315 boost::shared_ptr<MatrixDouble> mat_vals)
319 boost::shared_ptr<MatrixDouble> mat_vals,
322 betaCoeff(beta_fun) {}
343 template <
int SPACE_DIM, IntegrationType I,
typename OpBase>
346 template <
int SPACE_DIM,
typename OpBase>
349 boost::shared_ptr<MatrixDouble> mat_vals)
353 boost::shared_ptr<MatrixDouble> mat_vals,
356 betaCoeff(beta_fun) {}
377 template <
int SPACE_DIM, IntegrationType I,
typename OpBase>
386 template <
int FIELD_DIM, IntegrationType I,
typename OpBase>
394 template <
typename OpBase>
398 ScalarFun source_fun = [](
double,
double,
double) constexpr {
return 1; },
399 boost::shared_ptr<Range> ents_ptr =
nullptr)
401 sourceFun(source_fun) {}
409 template <
typename OpBase>
413 ScalarFun source_fun = [](
double,
double,
double) constexpr {
return 1; },
414 boost::shared_ptr<Range> ents_ptr =
nullptr)
416 sourceFun(source_fun) {}
434 template <
int SPACE_DIM, IntegrationType I,
typename OpBase>
437 template <
int SPACE_DIM,
typename OpBase>
441 const std::string
field_name, boost::shared_ptr<MatrixDouble> u_ptr,
442 ScalarFun beta_coeff = [](
double,
double,
double) constexpr {
return 1; },
443 boost::shared_ptr<Range> ents_ptr =
nullptr)
445 betaCoeff(beta_coeff) {}
448 boost::shared_ptr<MatrixDouble>
uPtr;
457 template <
int SPACE_DIM,
typename OpBase>
460 const std::string
field_name, boost::shared_ptr<MatrixDouble> u_ptr,
461 boost::shared_ptr<MatrixDouble> y_grad_ptr,
462 ConstantFun source_fun = []() constexpr {
return 1; })
464 yGradPtr(y_grad_ptr), alphaConstant(source_fun) {}
467 boost::shared_ptr<MatrixDouble>
uPtr;
473 template <
int FIELD_DIM,
int SPACE_DIM,
typename OpBase>
477 const std::string
field_name, boost::shared_ptr<MatrixDouble> u_ptr,
478 boost::shared_ptr<MatrixDouble> y_grad_ptr,
479 ConstantFun source_fun = []() constexpr {
return 1; })
481 yGradPtr(y_grad_ptr), alphaConstant(source_fun) {}
484 boost::shared_ptr<MatrixDouble>
uPtr;
498 template <
typename EleOp>
499 template <AssemblyType A>
500 template <IntegrationType I>
523 template <
int BASE_DIM,
int S = 1>
528 template <
int BASE_DIM,
int S = 1>
529 using OpBaseTimesScalarField = OpBaseTimesScalar<BASE_DIM, S>;
539 template <
int BASE_DIM,
int FIELD_DIM,
int S>
540 using OpBaseTimesVector =
556 template <
int BASE_DIM,
int FIELD_DIM,
int SPACE_DIM,
int S = 1>
570 template <
int BASE_DIM,
int FIELD_DIM,
int SPACE_DIM,
int S = 1>
571 using OpGradTimesSymTensor =
583 using OpMixDivTimesU =
591 template <
int SPACE_DIM>
599 template <
int SPACE_DIM>
600 using OpMixVecTimesDivLambda =
613 template <
int SPACE_DIM>
614 using OpNormalMixVecTimesScalar =
627 template <
int SPACE_DIM>
628 using OpNormalMixVecTimesVectorField =
643 template <
int BASE_DIM,
int FIELD_DIM,
int SPACE_DIM>
648 template <
typename OpBase>
655 const double vol = OpBase::getMeasure();
657 auto t_w = OpBase::getFTensor0IntegrationWeight();
661 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
666 t_w * vol * sourceFun(t_coords(0), t_coords(1), t_coords(2));
681 template <
int FIELD_DIM,
typename OpBase>
689 const double vol = OpBase::getMeasure();
691 auto t_w = OpBase::getFTensor0IntegrationWeight();
695 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
699 auto t_source = sourceFun(t_coords(0), t_coords(1), t_coords(2));
701 const double alpha = t_w * vol;
703 auto t_nf = getFTensor1FromArray<FIELD_DIM, FIELD_DIM>(
OpBase::locF);
706 t_nf(
i) += alpha * t_row_base * t_source(
i);
718 template <
int FIELD_DIM,
typename OpBase>
725 const size_t nb_base_functions = row_data.
getN().size2() / 3;
727 const double vol = OpBase::getMeasure();
729 auto t_w = OpBase::getFTensor0IntegrationWeight();
733 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
737 auto t_source = sourceFun(t_coords(0), t_coords(1), t_coords(2));
739 const double alpha = t_w * vol;
746 for (; rr < nb_base_functions; ++rr)
754 template <
int S,
typename OpBase>
760 const double vol = OpBase::getMeasure();
762 auto t_w = OpBase::getFTensor0IntegrationWeight();
766 auto t_vec = getFTensor0FromVec<S>(*sourceVec);
771 "Wrong number of integration points %d != %d",
777 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
782 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
798 template <
int FIELD_DIM,
int S,
typename OpBase>
804 const double vol = OpBase::getMeasure();
806 auto t_w = OpBase::getFTensor0IntegrationWeight();
810 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
812 auto t_vec = getFTensor1FromMat<FIELD_DIM, S>(*sourceVec);
817 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
819 auto t_nf = OpBase::template getNf<FIELD_DIM>();
823 t_nf(
i) += alpha * t_row_base * t_vec(
i);
836 template <
int FIELD_DIM,
int S,
typename OpBase>
841 const size_t nb_base_functions = row_data.
getN().size2() / 3;
843 const double vol = OpBase::getMeasure();
845 auto t_w = OpBase::getFTensor0IntegrationWeight();
849 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
851 auto t_vec = getFTensor1FromMat<FIELD_DIM, S>(*sourceVec);
856 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
863 for (; rr < nb_base_functions; ++rr)
872 template <
int SPACE_DIM,
int S,
typename OpBase>
879 const double vol = OpBase::getMeasure();
881 auto t_w = OpBase::getFTensor0IntegrationWeight();
885 auto t_val_grad = getFTensor1FromMat<SPACE_DIM, S>(*(matVals));
887 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
890 const double beta = vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
892 const double alpha = t_w * beta;
911 template <
int SPACE_DIM,
int S,
typename OpBase>
918 const double vol = OpBase::getMeasure();
920 auto t_w = OpBase::getFTensor0IntegrationWeight();
924 auto t_val_grad = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(matVals));
926 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
931 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
933 auto t_nf = OpBase::template getNf<SPACE_DIM>();
938 t_nf(
i) += alpha * (t_row_grad(
j) * t_val_grad(
i,
j));
953 template <
int SPACE_DIM,
int S,
typename OpBase>
959 const double vol = OpBase::getMeasure();
961 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
963 auto t_w = OpBase::getFTensor0IntegrationWeight();
967 auto t_val_mat = getFTensor2SymmetricFromMat<SPACE_DIM, S>(*(matVals));
972 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
974 auto t_nf = OpBase::template getNf<SPACE_DIM>();
979 t_nf(
j) += alpha * (t_row_grad(
i) * t_val_mat(
i,
j));
1000 const size_t nb_base_functions = row_data.
getN().size2() / 3;
1001 auto t_w = this->getFTensor0IntegrationWeight();
1003 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1006 auto t_u = getFTensor1FromMat<FIELD_DIM>(*(matVals));
1010 const double alpha = this->getMeasure() * t_w *
1011 betaConst(t_coords(0), t_coords(1), t_coords(2));
1012 auto t_nf = OpBase::template getNf<FIELD_DIM>();
1015 for (; bb != this->nbRows /
FIELD_DIM; ++bb) {
1016 const double t_div_base = t_diff_base(
j,
j);
1017 t_nf(
i) += alpha * t_div_base * t_u(
i);
1019 t_nf(
i) += t_base(0) * (alpha / t_coords(0)) * t_u(
i);
1025 for (; bb < nb_base_functions; ++bb) {
1038 template <
int SPACE_DIM,
typename OpBase, CoordinateTypes CoordSys>
1044 const size_t nb_base_functions = row_data.
getN().size2() / 3;
1045 auto t_w = this->getFTensor0IntegrationWeight();
1047 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1053 const double alpha = this->getMeasure() * t_w *
1054 betaConst(t_coords(0), t_coords(1), t_coords(2));
1058 for (; bb != this->nbRows; ++bb) {
1059 const double t_div_base = t_diff_base(
j,
j);
1063 for (; bb < nb_base_functions; ++bb)
1094 template <
int FIELD_DIM,
typename OpBase, CoordinateTypes CoordSys>
1102 const double vol = OpBase::getMeasure();
1104 auto t_w = OpBase::getFTensor0IntegrationWeight();
1106 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1114 const double alpha =
1115 t_w * vol * betaCoeff(t_coords(0), t_coords(1), t_coords(2));
1116 auto t_nf = getFTensor1FromArray<FIELD_DIM, FIELD_DIM>(
OpBase::locF);
1120 t_nf(
i) += alpha * t_row_grad(
i) * t_vec;
1133 template <
int SPACE_DIM,
typename OpBase>
1138 const size_t nb_base_functions = row_data.
getN().size2() / 3;
1139 auto t_w = this->getFTensor0IntegrationWeight();
1140 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1142 auto t_grad = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(*(matVals));
1146 const double alpha = this->getMeasure() * t_w;
1147 auto t_nf = OpBase::template getNf<SPACE_DIM>();
1150 for (; bb != this->nbRows /
SPACE_DIM; ++bb) {
1151 t_nf(
i) += alpha * betaCoeff(t_coords(0), t_coords(1), t_coords(2)) *
1152 t_base(
j) * t_grad(
i,
j);
1156 for (; bb < nb_base_functions; ++bb)
1167 template <
int SPACE_DIM,
typename OpBase>
1172 const size_t nb_base_functions = row_data.
getN().size2();
1173 auto t_w = this->getFTensor0IntegrationWeight();
1174 auto t_coords = this->getFTensor1CoordsAtGaussPts();
1176 auto t_div = getFTensor1FromMat<SPACE_DIM>(*(matVals));
1179 const double alpha = this->getMeasure() * t_w;
1180 auto t_nf = OpBase::template getNf<SPACE_DIM>();
1183 for (; bb != this->nbRows /
SPACE_DIM; ++bb) {
1184 t_nf(
i) += alpha * t_base *
1185 betaCoeff(t_coords(0), t_coords(1), t_coords(2)) * t_div(
i);
1189 for (; bb < nb_base_functions; ++bb)
1200 template <
typename OpBase>
1205 const size_t nb_base_functions = row_data.
getN().size2() / 3;
1208 auto t_w = OpBase::getFTensor0IntegrationWeight();
1212 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1214 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
1216 if (this->getNumeredEntFiniteElementPtr()->getEntType() == MBTRI)
1221 const double alpha =
1222 t_w * sourceFun(t_coords(0), t_coords(1), t_coords(2)) /
a;
1229 for (; rr < nb_base_functions; ++rr)
1238 template <
typename OpBase>
1243 const size_t nb_base_functions = row_data.
getN().size2() / 3;
1247 auto t_w = OpBase::getFTensor0IntegrationWeight();
1251 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1253 auto t_tangent = OpBase::getFTensor1TangentAtGaussPts();
1257 const double alpha = t_w * sourceFun(t_coords(0), t_coords(1), t_coords(2));
1267 for (; rr < nb_base_functions; ++rr)
1276 template <
int SPACE_DIM,
typename OpBase>
1284 const size_t nb_base_functions = row_data.
getN().size2() / 3;
1287 auto t_w = OpBase::getFTensor0IntegrationWeight();
1291 auto t_coords = OpBase::getFTensor1CoordsAtGaussPts();
1293 auto t_normal = OpBase::getFTensor1NormalsAtGaussPts();
1295 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
1297 auto a = OpBase::getMeasure();
1300 auto l2 = std::sqrt(t_normal(
i) * t_normal(
i));
1301 const double alpha =
1302 t_w * betaCoeff(t_coords(0), t_coords(1), t_coords(2)) * (
a / l2);
1304 auto t_nf = OpBase::template getNf<SPACE_DIM>();
1308 t_nf(
J) += alpha * (t_row_base(
i) * t_normal(
i)) * t_u(
J);
1312 for (; rr < nb_base_functions; ++rr)
1322 template <
int SPACE_DIM,
typename OpBase>
1328 auto t_w = this->getFTensor0IntegrationWeight();
1331 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
1332 auto t_grad_y = getFTensor1FromMat<SPACE_DIM>(*yGradPtr);
1335 const double alpha_constant = alphaConstant();
1339 const double vol = OpBase::getMeasure();
1340 const double c = (t_grad_y(
i) * t_u(
i)) * (t_w * vol * alpha_constant);
1359 template <
int FIELD_DIM,
int SPACE_DIM,
typename OpBase>
1365 auto t_w = this->getFTensor0IntegrationWeight();
1368 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
1369 auto t_grad_y = getFTensor2FromMat<FIELD_DIM, SPACE_DIM>(*yGradPtr);
1373 const double alpha_constant = alphaConstant();
1377 const double vol = OpBase::getMeasure();
1380 t_c(
J) = (t_grad_y(
J,
i) * t_u(
i)) * (t_w * vol * alpha_constant);
1382 auto t_nf = getFTensor1FromArray<FIELD_DIM, FIELD_DIM>(
OpBase::locF);
1385 t_nf(
J) += t_base * t_c(
J);
1402 #endif // __LINEAR_FORMS_INTEGRATORS_HPP__