v0.14.0
Public Member Functions | Public Attributes | List of all members
ConvectiveMassElement::OpMassJacobian Struct Reference

#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

BlockDatadAta
 
CommonDatacommonData
 
int tAg
 
bool jAcobian
 
boollInear
 
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< doubleactive
 

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 154 of file ConvectiveMassElement.cpp.

159  field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
160  dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
161  lInear(commonData.lInear), fieldDisp(false), methodsOp(methods_op) {}

Member Function Documentation

◆ doWork()

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

Definition at line 163 of file ConvectiveMassElement.cpp.

164  {
166 
167  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
168  dAta.tEts.end()) {
170  }
171 
172  // do it only once, no need to repeat this for edges,faces or tets
173  if (row_type != MBVERTEX)
175 
176  int nb_dofs = row_data.getIndices().size();
177  if (nb_dofs == 0)
179 
180  {
181 
182  if (a.size() != 3) {
183  a.resize(3, false);
184  dot_W.resize(3, false);
185  a_res.resize(3, false);
186  g.resize(3, 3, false);
187  G.resize(3, 3, false);
188  h.resize(3, 3, false);
189  H.resize(3, 3, false);
190  invH.resize(3, 3, false);
191  F.resize(3, 3, false);
192  }
193 
194  std::fill(dot_W.begin(), dot_W.end(), 0);
195  std::fill(H.data().begin(), H.data().end(), 0);
196  std::fill(invH.data().begin(), invH.data().end(), 0);
197  for (int ii = 0; ii != 3; ii++) {
198  H(ii, ii) = 1;
199  invH(ii, ii) = 1;
200  }
201 
202  int nb_gauss_pts = row_data.getN().size1();
203  commonData.valMass.resize(nb_gauss_pts);
204  commonData.jacMassRowPtr.resize(nb_gauss_pts);
205  commonData.jacMass.resize(nb_gauss_pts);
206 
207  const std::vector<VectorDouble> &dot_spacial_vel =
209 
210  const std::vector<MatrixDouble> &spatial_positions_grad =
212 
213  const std::vector<MatrixDouble> &spatial_velocities_grad =
215 
216  const std::vector<VectorDouble> &meshpos_vel =
218 
219  const std::vector<MatrixDouble> &mesh_positions_gradient =
221 
222  int nb_active_vars = 0;
223  for (int gg = 0; gg < nb_gauss_pts; gg++) {
224 
225  if (gg == 0) {
226 
227  trace_on(tAg);
228 
229  for (int nn1 = 0; nn1 < 3; nn1++) { // 0
230  // commonData.dataAtGaussPts["DOT_"+commonData.spatialVelocities]
231  a[nn1] <<= dot_spacial_vel[gg][nn1];
232  nb_active_vars++;
233  }
234  for (int nn1 = 0; nn1 < 3; nn1++) { // 3
235  for (int nn2 = 0; nn2 < 3; nn2++) {
236  // commonData.gradAtGaussPts[commonData.spatialPositions][gg]
237  h(nn1, nn2) <<= spatial_positions_grad[gg](nn1, nn2);
238  if (fieldDisp) {
239  if (nn1 == nn2) {
240  h(nn1, nn2) += 1;
241  }
242  }
243  nb_active_vars++;
244  }
245  }
247  .size() > 0) {
248  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
249  for (int nn2 = 0; nn2 < 3; nn2++) {
250  // commonData.gradAtGaussPts[commonData.spatialVelocities]
251  g(nn1, nn2) <<= spatial_velocities_grad[gg](nn1, nn2);
252  nb_active_vars++;
253  }
254  }
255  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
256  // commonData.dataAtGaussPts["DOT_"+commonData.meshPositions]
257  dot_W(nn1) <<= meshpos_vel[gg][nn1];
258  nb_active_vars++;
259  }
260  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
261  for (int nn2 = 0; nn2 < 3; nn2++) {
262  // commonData.gradAtGaussPts[commonData.meshPositions][gg]
263  H(nn1, nn2) <<= mesh_positions_gradient[gg](nn1, nn2);
264  nb_active_vars++;
265  }
266  }
267  }
268 
269  auto a0 = dAta.a0;
271 
272  auto t_a_res =
274  auto t_a = FTensor::Tensor1<adouble *, 3>{&a[0], &a[1], &a[2]};
275  auto t_a0 = FTensor::Tensor1<double *, 3>{&a0[0], &a0[1], &a0[2]};
276  auto t_dotW =
280  auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
283 
284  const double rho0 = dAta.rho0;
285 
287  CHKERR invertTensor3by3(H, detH, invH);
288 
289  t_G(i, j) = t_g(i, k) * t_invH(k, j);
290  t_a_res(i) = t_a(i) - t_a0(i) + t_G(i, j) * t_dotW(j);
291 
292  // FIXME: there is error somewhere for nonlinear case
293  // test dam example with -is_linear 0
294  if (!lInear) {
295 
296  t_F(i, j) = t_h(i, k) * t_invH(k, j);
297  t_a_res(i) *= rho0 * detH;
298  t_a_res(i) *= determinantTensor3by3(t_F);
299 
300  } else {
301 
302  t_a_res(i) *= rho0 * detH;
303  }
304 
305  // dependant
306  VectorDouble &res = commonData.valMass[gg];
307  res.resize(3);
308  for (int rr = 0; rr < 3; rr++) {
309  a_res[rr] >>= res[rr];
310  }
311 
312  trace_off();
313  }
314 
315  active.resize(nb_active_vars);
316  int aa = 0;
317  for (int nn1 = 0; nn1 < 3; nn1++) { // 0
318  active[aa++] = dot_spacial_vel[gg][nn1];
319  }
320  for (int nn1 = 0; nn1 < 3; nn1++) { // 3
321  for (int nn2 = 0; nn2 < 3; nn2++) {
322  if (fieldDisp && nn1 == nn2) {
323  active[aa++] = spatial_positions_grad[gg](nn1, nn2) + 1;
324  } else {
325  active[aa++] = spatial_positions_grad[gg](nn1, nn2);
326  }
327  }
328  }
329  if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() >
330  0) {
331  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9=12
332  for (int nn2 = 0; nn2 < 3; nn2++) {
333  active[aa++] = spatial_velocities_grad[gg](nn1, nn2);
334  }
335  }
336  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9=21
337  active[aa++] = meshpos_vel[gg][nn1];
338  }
339  for (int nn1 = 0; nn1 < 3; nn1++) { // 3+9+9+3=24
340  for (int nn2 = 0; nn2 < 3; nn2++) {
341  active[aa++] = mesh_positions_gradient[gg](nn1, nn2);
342  }
343  }
344  }
345 
346  if (!jAcobian) {
347  VectorDouble &res = commonData.valMass[gg];
348  if (gg > 0) {
349  res.resize(3);
350  int r;
351  r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
352  if (r != 3) { // function is locally analytic
353  SETERRQ1(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
354  "ADOL-C function evaluation with error r = %d", r);
355  }
356  }
357  double val = getVolume() * getGaussPts()(3, gg);
358  res *= val;
359  // cout << "my res " << res << endl;
360  } else {
361  commonData.jacMassRowPtr[gg].resize(3);
362  commonData.jacMass[gg].resize(3, nb_active_vars);
363  for (int nn1 = 0; nn1 < 3; nn1++) {
364  (commonData.jacMassRowPtr[gg])[nn1] =
365  &(commonData.jacMass[gg](nn1, 0));
366  }
367  int r;
368  r = jacobian(tAg, 3, nb_active_vars, &active[0],
369  &(commonData.jacMassRowPtr[gg])[0]);
370  if (r != 3) {
371  SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
372  "ADOL-C function evaluation with error");
373  }
374  double val = getVolume() * getGaussPts()(3, gg);
375  commonData.jacMass[gg] *= val;
376  }
377  }
378  }
379 
381 }

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:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
ConvectiveMassElement::OpMassJacobian::fieldDisp
bool fieldDisp
Definition: ConvectiveMassElement.hpp:186
ConvectiveMassElement::OpMassJacobian::active
std::vector< double > active
Definition: ConvectiveMassElement.hpp:201
ConvectiveMassElement::OpMassJacobian::j
FTensor::Index< 'j', 3 > j
Definition: ConvectiveMassElement.hpp:196
ConvectiveMassElement::BlockData::tEts
Range tEts
elements in block set
Definition: ConvectiveMassElement.hpp:119
ConvectiveMassElement::OpMassJacobian::dot_W
VectorBoundedArray< adouble, 3 > dot_W
Definition: ConvectiveMassElement.hpp:199
ConvectiveMassElement::CommonData::spatialVelocities
string spatialVelocities
Definition: ConvectiveMassElement.hpp:137
ConvectiveMassElement::CommonData::dataAtGaussPts
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
Definition: ConvectiveMassElement.hpp:133
sdf.r
int r
Definition: sdf.py:8
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:1588
ConvectiveMassElement::CommonData::jacMassRowPtr
std::vector< std::vector< double * > > jacMassRowPtr
Definition: ConvectiveMassElement.hpp:142
ConvectiveMassElement::CommonData::gradAtGaussPts
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Definition: ConvectiveMassElement.hpp:134
FTensor::Number< 0 >
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
ConvectiveMassElement::OpMassJacobian::G
MatrixBoundedArray< adouble, 9 > G
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::a
VectorBoundedArray< adouble, 3 > a
Definition: ConvectiveMassElement.hpp:199
ConvectiveMassElement::OpMassJacobian::i
FTensor::Index< 'i', 3 > i
Definition: ConvectiveMassElement.hpp:195
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1214
ConvectiveMassElement::OpMassJacobian::H
MatrixBoundedArray< adouble, 9 > H
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::k
FTensor::Index< 'k', 3 > k
Definition: ConvectiveMassElement.hpp:197
ConvectiveMassElement::OpMassJacobian::jAcobian
bool jAcobian
Definition: ConvectiveMassElement.hpp:184
ConvectiveMassElement::CommonData::jacMass
std::vector< MatrixDouble > jacMass
Definition: ConvectiveMassElement.hpp:143
ConvectiveMassElement::OpMassJacobian::h
MatrixBoundedArray< adouble, 9 > h
Definition: ConvectiveMassElement.hpp:200
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
ConvectiveMassElement::OpMassJacobian::methodsOp
boost::ptr_vector< MethodForForceScaling > & methodsOp
Definition: ConvectiveMassElement.hpp:188
ConvectiveMassElement::CommonData::spatialPositions
string spatialPositions
Definition: ConvectiveMassElement.hpp:135
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
ConvectiveMassElement::CommonData::lInear
bool lInear
Definition: ConvectiveMassElement.hpp:129
ConvectiveMassElement::OpMassJacobian::g
MatrixBoundedArray< adouble, 9 > g
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::BlockData::rho0
double rho0
reference density
Definition: ConvectiveMassElement.hpp:117
ConvectiveMassElement::OpMassJacobian::commonData
CommonData & commonData
Definition: ConvectiveMassElement.hpp:182
ConvectiveMassElement::OpMassJacobian::invH
MatrixBoundedArray< adouble, 9 > invH
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::a_res
VectorBoundedArray< adouble, 3 > a_res
Definition: ConvectiveMassElement.hpp:199
MoFEM::getFTensor2FromArray3by3
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr, const size_t cc=0)
Definition: Templates.hpp:1342
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
adouble
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:1318
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
ConvectiveMassElement::OpMassJacobian::F
MatrixBoundedArray< adouble, 9 > F
Definition: ConvectiveMassElement.hpp:200
ConvectiveMassElement::OpMassJacobian::lInear
bool & lInear
Definition: ConvectiveMassElement.hpp:185
ConvectiveMassElement::CommonData::valMass
std::vector< VectorDouble > valMass
Definition: ConvectiveMassElement.hpp:141
ConvectiveMassElement::OpMassJacobian::dAta
BlockData & dAta
Definition: ConvectiveMassElement.hpp:181
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
a0
constexpr double a0
Definition: hcurl_check_approx_in_2d.cpp:37
ConvectiveMassElement::CommonData::meshPositions
string meshPositions
Definition: ConvectiveMassElement.hpp:136
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MethodForForceScaling::applyScale
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Definition: MethodForForceScaling.hpp:21
ConvectiveMassElement::OpMassJacobian::tAg
int tAg
Definition: ConvectiveMassElement.hpp:183
ConvectiveMassElement::BlockData::a0
VectorDouble a0
constant acceleration
Definition: ConvectiveMassElement.hpp:118