13 boost::shared_ptr<MatrixDouble> base_mass_ptr,
14 boost::shared_ptr<EntitiesFieldData> data_l2,
17 verbosity(verb), severityLevel(sev), baseMassPtr(base_mass_ptr),
27 "Space should be set to L2");
30 auto fe_type = getFEType();
32 auto calculate_base = [&]() {
34 auto fe_ptr = getPtrFE();
36 dataL2->dataOnEntities[fe_type].clear();
37 dataL2->dataOnEntities[MBVERTEX].push_back(
41 auto &vertex_data = dataL2->dataOnEntities[MBVERTEX][0];
42 vertex_data.getNSharedPtr(
NOBASE) =
43 fe_ptr->getEntData(
H1, MBVERTEX, 0).getNSharedPtr(
NOBASE);
44 vertex_data.getDiffNSharedPtr(
NOBASE) =
45 fe_ptr->getEntData(
H1, MBVERTEX, 0).getDiffNSharedPtr(
NOBASE);
47 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
48 ent_data.getSense() = 1;
49 ent_data.getBase() = base;
50 ent_data.getOrder() = std::max(0, fe_ptr->getMaxDataOrder() - 1);
52 CHKERR fe_ptr->getElementPolynomialBase()->getValue(
53 getGaussPts(), boost::make_shared<EntPolynomialBaseCtx>(
61 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
62 auto &base_funcions = ent_data.getN(base);
63 const auto nb = base_funcions.size2();
67 const auto nb_integration_pts = getGaussPts().size2();
72 "Mass matrix is null pointer");
75 auto &nN = *baseMassPtr;
76 nN.resize(nb, nb,
false);
79 auto t_w = getFTensor0IntegrationWeight();
81 auto t_row_base = ent_data.getFTensor0N();
83 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
85 const double alpha = t_w;
87 for (
int rr = 0; rr != nb; ++rr) {
90 auto a_mat_ptr = &nN(rr, 0);
92 auto t_col_base = ent_data.getFTensor0N(gg, 0);
94 for (
int cc = 0; cc <= rr; ++cc) {
96 *a_mat_ptr += alpha * (t_row_base * t_col_base);
113 boost::shared_ptr<EntitiesFieldData> data_l2,
114 boost::shared_ptr<MatrixDouble> inv_jac_ptr)
124 return applyTransform<2, 2, 2, 2>(diff_n);
127 const auto fe_type = getFEType();
128 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
129 CHKERR apply_transform(ent_data.getDiffN());
135 int derivative, boost::shared_ptr<MatrixDouble> base_mass_ptr,
136 boost::shared_ptr<EntitiesFieldData> data_l2,
139 calcBaseDerivative(derivative) {
141 doEntities[MBVERTEX] =
false;
144template <
int BASE_DIM,
int SPACE_DIM>
150 const int nb_gauss_pts = data.
getN(base).size1();
151 const int nb_approx_bases = data.
getN(base).size2() /
BASE_DIM;
152 const int nb_derivatives =
155 const int nb_prj_bases = ent_data.
getN().size2();
159 if (!n_diff_shared_ptr)
160 n_diff_shared_ptr = boost::make_shared<MatrixDouble>();
162 auto &nex_diff_base = *(n_diff_shared_ptr);
164 nex_diff_base.resize(nb_gauss_pts, nb_approx_bases * next_nb_derivatives,
166 nex_diff_base.clear();
169 auto next_diffs_ptr = &*nex_diff_base.data().begin();
170 auto t_next_diff = getFTensor1FromPtr<SPACE_DIM>(next_diffs_ptr);
172 for (
int gg = 0; gg != nb_gauss_pts; ++gg) {
174 auto ptr = &*nF.data().begin();
176 for (
auto r = 0; r != nb_approx_bases * nb_derivatives; ++r) {
179 for (
int rr = 0; rr != nb_prj_bases; ++rr) {
180 t_next_diff(
i) += l2_diff_base(
i) * (*ptr);
192template <
int BASE_DIM>
198 auto &approx_base = data.
getN(base);
199 const auto nb_approx_bases = approx_base.size2() /
BASE_DIM;
201 if (nb_approx_bases) {
203 const auto fe_type = getFEType();
204 const auto nb_integration_pts = approx_base.size1();
206 const auto space_dim =
210 int nb_derivatives =
BASE_DIM * pow(space_dim, calcBaseDerivative - 1);
212 auto &ent_data = dataL2->dataOnEntities[fe_type][0];
213 const int nb_prj_bases = ent_data.getN().size2();
218 "Mass matrix is null pointer");
220 auto &nN = *baseMassPtr;
223 if (diff_approx_base.size2() != nb_approx_bases * nb_derivatives) {
225 "Number of derivatives and basses do not match");
227 if (ent_data.getN().size1() != nb_integration_pts) {
229 "Number of integration points is not consistent");
231 if (nN.size2() != nb_prj_bases) {
233 "Number of base functions and size of mass matrix does not math");
237 nF.resize(nb_approx_bases * nb_derivatives, nb_prj_bases,
false);
240 auto t_w = getFTensor0IntegrationWeight();
243 auto diff_base_ptr = &*diff_approx_base.data().begin();
245 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
247 const double alpha = t_w;
249 for (
int r = 0; r != nb_approx_bases * nb_derivatives; ++r) {
252 auto t_base = ent_data.getFTensor0N(base, gg, 0);
253 for (
int rr = 0; rr != nb_prj_bases; ++rr) {
254 nF(r, rr) += alpha * (t_base * (*diff_base_ptr));
264 for (
auto r = 0; r != nb_approx_bases * nb_derivatives; ++r) {
265 ublas::matrix_row<MatrixDouble> mc(nF, r);
270 CHKERR setBaseImpl<BASE_DIM, 3>(data, ent_data);
271 else if (space_dim == 2)
272 CHKERR setBaseImpl<BASE_DIM, 2>(data, ent_data);
277 "Space dim can be only 1,2,3 but is %d", space_dim);
286 return doWorkImpl<1>(side, type, data);
292 return doWorkImpl<3>(side, type, data);
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
FieldApproximationBase
approximation base
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
FieldSpace
approximation spaces
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
SeverityLevel
Severity levels.
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives direvatie)
structure to get information form mofem into EntitiesFieldData
OpBaseDerivativesBase(boost::shared_ptr< MatrixDouble > base_mass_ptr, boost::shared_ptr< EntitiesFieldData > data_l2, const FieldApproximationBase b, const FieldSpace s, int verb=QUIET, Sev sev=Sev::verbose)
Transform local reference derivatives of shape functions to global derivatives.