v0.13.1
Public Member Functions | Public Attributes | List of all members
ConvectiveMassElement::OpMassJacobian Struct Reference
User modules » Mass Element

#include <users_modules/basic_finite_elements/src/ConvectiveMassElement.hpp>

Inheritance diagram for ConvectiveMassElement::OpMassJacobian:
[legend]
Collaboration diagram for ConvectiveMassElement::OpMassJacobian:
[legend]

Public Member Functions

 OpMassJacobian (const std::string field_name, BlockData &data, CommonData &common_data, boost::ptr_vector< MethodForForceScaling > &methods_op, int tag, bool linear=false)
 
MoFEMErrorCode doWork (int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)
 

Public Attributes

BlockData & dAta
 
CommonData & commonData
 
int tAg
 
bool jAcobian
 
bool & lInear
 
bool fieldDisp
 
boost::ptr_vector< MethodForForceScaling > & methodsOp
 
FTensor::Index< 'i', 3 > i
 
FTensor::Index< 'j', 3 > j
 
FTensor::Index< 'k', 3 > k
 
VectorBoundedArray< adouble, 3 > a
 
VectorBoundedArray< adouble, 3 > dot_W
 
VectorBoundedArray< adouble, 3 > dp_dt
 
VectorBoundedArray< adouble, 3 > a_res
 
MatrixBoundedArray< adouble, 9 > h
 
MatrixBoundedArray< adouble, 9 > H
 
MatrixBoundedArray< adouble, 9 > invH
 
MatrixBoundedArray< adouble, 9 > F
 
MatrixBoundedArray< adouble, 9 > g
 
MatrixBoundedArray< adouble, 9 > G
 
std::vector< double > active
 

Detailed Description

Definition at line 177 of file ConvectiveMassElement.hpp.

Constructor & Destructor Documentation

◆ OpMassJacobian()

ConvectiveMassElement::OpMassJacobian::OpMassJacobian ( const std::string  field_name,
BlockData &  data,
CommonData &  common_data,
boost::ptr_vector< MethodForForceScaling > &  methods_op,
int  tag,
bool  linear = false 
)

Definition at line 156 of file ConvectiveMassElement.cpp.

160 : VolumeElementForcesAndSourcesCore::UserDataOperator(
161 field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
162 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
163 lInear(commonData.lInear), fieldDisp(false), methodsOp(methods_op) {}
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:17
ConvectiveMassElement::OpMassJacobian::methodsOp
boost::ptr_vector< MethodForForceScaling > & methodsOp
Definition: ConvectiveMassElement.hpp:188

Member Function Documentation

◆ doWork()

MoFEMErrorCode ConvectiveMassElement::OpMassJacobian::doWork ( int  row_side,
EntityType  row_type,
EntitiesFieldData::EntData &  row_data 
)

Definition at line 165 of file ConvectiveMassElement.cpp.

167 {
168 MoFEMFunctionBeginHot;
169
170 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
171 dAta.tEts.end()) {
172 MoFEMFunctionReturnHot(0);
173 }
174
175 // do it only once, no need to repeat this for edges,faces or tets
176 if (row_type != MBVERTEX)
177 MoFEMFunctionReturnHot(0);
178
179 int nb_dofs = row_data.getIndices().size();
180 if (nb_dofs == 0)
181 MoFEMFunctionReturnHot(0);
182
183 {
184
185 if (a.size() != 3) {
186 a.resize(3, false);
187 dot_W.resize(3, false);
188 a_res.resize(3, false);
189 g.resize(3, 3, false);
190 G.resize(3, 3, false);
191 h.resize(3, 3, false);
192 H.resize(3, 3, false);
193 invH.resize(3, 3, false);
194 F.resize(3, 3, false);
195 }
196
197 dot_W.clear();
198 H.clear();
199 invH.clear();
200 for (int dd = 0; dd < 3; dd++) {
201 H(dd, dd) = 1;
202 invH(dd, dd) = 1;
203 }
204
205 int nb_gauss_pts = row_data.getN().size1();
206 commonData.valMass.resize(nb_gauss_pts);
207 commonData.jacMassRowPtr.resize(nb_gauss_pts);
208 commonData.jacMass.resize(nb_gauss_pts);
209
210 const std::vector<VectorDouble> &dot_spacial_vel =
211 commonData.dataAtGaussPts["DOT_" + commonData.spatialVelocities];
212
213 const std::vector<MatrixDouble> &spatial_positions_grad =
214 commonData.gradAtGaussPts[commonData.spatialPositions];
215
216 const std::vector<MatrixDouble> &spatial_velocities_grad =
217 commonData.gradAtGaussPts[commonData.spatialVelocities];
218
219 const std::vector<VectorDouble> &meshpos_vel =
220 commonData.dataAtGaussPts["DOT_" + commonData.meshPositions];
221
222 const std::vector<MatrixDouble> &mesh_positions_gradient =
223 commonData.gradAtGaussPts[commonData.meshPositions];
224
225 int nb_active_vars = 0;
226 for (int gg = 0; gg < nb_gauss_pts; gg++) {
227
228 if (gg == 0) {
229
230 trace_on(tAg);
231
232 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
233 // commonData.dataAtGaussPts["DOT_"+commonData.spatialVelocities]
234 a[nn1] <<= dot_spacial_vel[gg][nn1];
235 nb_active_vars++;
236 }
237 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
238 for (int nn2 = 0; nn2 < 3; nn2++) {
239 // commonData.gradAtGaussPts[commonData.spatialPositions][gg]
240 h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
241 if (fieldDisp) {
242 if (nn1 == nn2) {
243 h(nn1, nn2) += 1;
244 }
245 }
246 nb_active_vars++;
247 }
248 }
249 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions]
250 .size() > 0) {
251 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
252 for (int nn2 = 0; nn2 < 3; nn2++) {
253 // commonData.gradAtGaussPts[commonData.spatialVelocities]
254 g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
255 nb_active_vars++;
256 }
257 }
258 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
259 // commonData.dataAtGaussPts["DOT_"+commonData.meshPositions]
260 dot_W(nn1) <<= meshpos_vel[gg][nn1];
261 nb_active_vars++;
262 }
263 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
264 for (int nn2 = 0; nn2 < 3; nn2++) {
265 // commonData.gradAtGaussPts[commonData.meshPositions][gg]
266 H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
267 nb_active_vars++;
268 }
269 }
270 }
271
272 auto a0 = dAta.a0;
273 CHKERR MethodForForceScaling::applyScale(getFEMethod(), methodsOp, a0);
274
275 auto t_a_res =
276 FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
277 auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
278 auto t_a0 = FTensor::Tensor1<double *, 3>{&a0[0], &a0[1], &a0[2]};
279 auto t_dotW =
280 FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
281 auto t_g = getFTensor2FromArray3by3(g, FTensor::Number<0>(), 0);
282 auto t_G = getFTensor2FromArray3by3(G, FTensor::Number<0>(), 0);
283 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
284 auto t_F = getFTensor2FromArray3by3(F, FTensor::Number<0>(), 0);
285 auto t_h = getFTensor2FromArray3by3(h, FTensor::Number<0>(), 0);
286
287 const double rho0 = dAta.rho0;
288
289 adouble detH = determinantTensor3by3(H);
290 CHKERR invertTensor3by3(H, detH, invH);
291
292 t_G(i, j) = t_g(i, k) * t_invH(k, j);
293 t_a_res(i) = t_a(i) - t_a0(i) + t_G(i, j) * t_dotW(j);
294
295 if (!lInear) {
296
297 t_F(i,j) = t_h(i,k)*t_invH(k,j);
298 t_a_res(i) *= rho0 * detH;
299 t_a_res(i) *= determinantTensor3by3(t_F);
300
301 } else {
302
303 t_a_res(i) *= rho0 * detH;
304
305 }
306
307 // dependant
308 VectorDouble &res = commonData.valMass[gg];
309 res.resize(3);
310 for (int rr = 0; rr < 3; rr++) {
311 a_res[rr] >>= res[rr];
312 }
313
314 trace_off();
315 }
316
317 active.resize(nb_active_vars);
318 int aa = 0;
319 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
320 active[aa++] = dot_spacial_vel[gg][nn1];
321 }
322 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
323 for (int nn2 = 0; nn2 < 3; nn2++) {
324 if (fieldDisp && nn1 == nn2) {
325 active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
326 } else {
327 active[aa++] = spatial_positions_grad[gg](nn1, nn2);
328 }
329 }
330 }
331 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() >
332 0) {
333 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
334 for (int nn2 = 0; nn2 < 3; nn2++) {
335 active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
336 }
337 }
338 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
339 active[aa++] = meshpos_vel[gg][nn1];
340 }
341 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
342 for (int nn2 = 0; nn2 < 3; nn2++) {
343 active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
344 }
345 }
346 }
347
348 if (!jAcobian) {
349 VectorDouble &res = commonData.valMass[gg];
350 if (gg > 0) {
351 res.resize(3);
352 int r;
353 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
354 if (r != 3) { // function is locally analytic
355 SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
356 "ADOL-C function evaluation with error r = %d", r);
357 }
358 }
359 double val = getVolume() * getGaussPts()(3, gg);
360 res *= val;
361 } else {
362 commonData.jacMassRowPtr[gg].resize(3);
363 commonData.jacMass[gg].resize(3, nb_active_vars);
364 for (int nn1 = 0; nn1 < 3; nn1++) {
365 (commonData.jacMassRowPtr[gg])[nn1] =
366 &(commonData.jacMass[gg](nn1, 0));
367 }
368 int r;
369 r = jacobian(tAg, 3, nb_active_vars, &active[0],
370 &(commonData.jacMassRowPtr[gg])[0]);
371 if (r != 3) {
372 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
373 "ADOL-C function evaluation with error");
374 }
375 double val = getVolume() * getGaussPts()(3, gg);
376 commonData.jacMass[gg] *= val;
377 }
378 }
379 }
380
381 MoFEMFunctionReturnHot(0);
382}
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
a0
constexpr double a0
Definition: hcurl_check_approx_in_2d.cpp:37
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
MoFEM::getFTensor2FromArray3by3
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
Definition: Templates.hpp:972
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1204
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1193
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:21
ConvectiveMassElement::BlockData::a0
VectorDouble a0
constant acceleration
Definition: ConvectiveMassElement.hpp:118
ConvectiveMassElement::CommonData::dataAtGaussPts
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
Definition: ConvectiveMassElement.hpp:133
ConvectiveMassElement::CommonData::gradAtGaussPts
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Definition: ConvectiveMassElement.hpp:134
ConvectiveMassElement::CommonData::jacMassRowPtr
std::vector< std::vector< double * > > jacMassRowPtr
Definition: ConvectiveMassElement.hpp:142
ConvectiveMassElement::OpMassJacobian::invH
MatrixBoundedArray< adouble, 9 > invH
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::G
MatrixBoundedArray< adouble, 9 > G
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::H
MatrixBoundedArray< adouble, 9 > H
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::g
MatrixBoundedArray< adouble, 9 > g
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::h
MatrixBoundedArray< adouble, 9 > h
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::a_res
VectorBoundedArray< adouble, 3 > a_res
Definition: ConvectiveMassElement.hpp:199
ConvectiveMassElement::OpMassJacobian::F
MatrixBoundedArray< adouble, 9 > F
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::dot_W
VectorBoundedArray< adouble, 3 > dot_W
Definition: ConvectiveMassElement.hpp:199
ConvectiveMassElement::OpMassJacobian::a
VectorBoundedArray< adouble, 3 > a
Definition: ConvectiveMassElement.hpp:199
MethodForForceScaling::applyScale
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Definition: MethodForForceScaling.hpp:21
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:1305
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201

Member Data Documentation

◆ a

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpMassJacobian::a

Definition at line 199 of file ConvectiveMassElement.hpp.

◆ a_res

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpMassJacobian::a_res

Definition at line 199 of file ConvectiveMassElement.hpp.

◆ active

std::vector<double> ConvectiveMassElement::OpMassJacobian::active

Definition at line 201 of file ConvectiveMassElement.hpp.

◆ commonData

CommonData& ConvectiveMassElement::OpMassJacobian::commonData

Definition at line 182 of file ConvectiveMassElement.hpp.

◆ dAta

BlockData& ConvectiveMassElement::OpMassJacobian::dAta

Definition at line 181 of file ConvectiveMassElement.hpp.

◆ dot_W

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpMassJacobian::dot_W

Definition at line 199 of file ConvectiveMassElement.hpp.

◆ dp_dt

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpMassJacobian::dp_dt

Definition at line 199 of file ConvectiveMassElement.hpp.

◆ F

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::F

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ fieldDisp

bool ConvectiveMassElement::OpMassJacobian::fieldDisp

Definition at line 186 of file ConvectiveMassElement.hpp.

◆ g

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::g

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ G

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::G

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ h

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::h

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ H

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::H

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ i

FTensor::Index<'i', 3> ConvectiveMassElement::OpMassJacobian::i

Definition at line 195 of file ConvectiveMassElement.hpp.

◆ invH

MatrixBoundedArray<adouble, 9> ConvectiveMassElement::OpMassJacobian::invH

Definition at line 200 of file ConvectiveMassElement.hpp.

◆ j

FTensor::Index<'j', 3> ConvectiveMassElement::OpMassJacobian::j

Definition at line 196 of file ConvectiveMassElement.hpp.

◆ jAcobian

bool ConvectiveMassElement::OpMassJacobian::jAcobian

Definition at line 184 of file ConvectiveMassElement.hpp.

◆ k

FTensor::Index<'k', 3> ConvectiveMassElement::OpMassJacobian::k

Definition at line 197 of file ConvectiveMassElement.hpp.

◆ lInear

bool& ConvectiveMassElement::OpMassJacobian::lInear

Definition at line 185 of file ConvectiveMassElement.hpp.

◆ methodsOp

boost::ptr_vector<MethodForForceScaling>& ConvectiveMassElement::OpMassJacobian::methodsOp

Definition at line 188 of file ConvectiveMassElement.hpp.

◆ tAg

int ConvectiveMassElement::OpMassJacobian::tAg

Definition at line 183 of file ConvectiveMassElement.hpp.

The documentation for this struct was generated from the following files:
Generated on Thu Aug 11 2022 23:36:35 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow