ConvectiveMassElement::OpVelocityJacobian Struct Reference

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

Inheritance diagram for ConvectiveMassElement::OpVelocityJacobian:

Collaboration diagram for ConvectiveMassElement::OpVelocityJacobian:

Public Member Functions
	OpVelocityJacobian (const std::string field_name, BlockData &data, CommonData &common_data, int tag, bool jacobian=true)
MoFEMErrorCode	doWork (int row_side, EntityType row_type, EntitiesFieldData::EntData &row_data)

Public Attributes
BlockData &	dAta
CommonData &	commonData
int	tAg
bool	jAcobian
bool	fieldDisp
VectorBoundedArray< adouble, 3 >	a_res
VectorBoundedArray< adouble, 3 >	v
VectorBoundedArray< adouble, 3 >	dot_w
VectorBoundedArray< adouble, 3 >	dot_W
VectorBoundedArray< adouble, 3 >	dot_u
MatrixBoundedArray< adouble, 9 >	h
MatrixBoundedArray< adouble, 9 >	H
MatrixBoundedArray< adouble, 9 >	invH
MatrixBoundedArray< adouble, 9 >	F
adouble	detH
std::vector< double >	active

Detailed Description

Definition at line 279 of file ConvectiveMassElement.hpp.

Constructor & Destructor Documentation

OpVelocityJacobian()

ConvectiveMassElement::OpVelocityJacobian::OpVelocityJacobian	(	const std::string	field_name,
		BlockData &	data,
		CommonData &	common_data,
		int	tag,
		bool	jacobian = true
	)

Definition at line 788 of file ConvectiveMassElement.cpp.

    : VolumeElementForcesAndSourcesCore::UserDataOperator(
          field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
      dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
      fieldDisp(false) {}

Member Function Documentation

doWork()

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

Definition at line 796 of file ConvectiveMassElement.cpp.

                                                                           {
  MoFEMFunctionBeginHot;
 
  if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
      dAta.tEts.end()) {
    MoFEMFunctionReturnHot(0);
  }
 
  // do it only once, no need to repeat this for edges,faces or tets
  if (row_type != MBVERTEX)
    MoFEMFunctionReturnHot(0);
 
  int nb_dofs = row_data.getIndices().size();
  if (nb_dofs == 0)
    MoFEMFunctionReturnHot(0);
 
  {
 
    v.resize(3);
    dot_w.resize(3);
    h.resize(3, 3);
    h.clear();
    F.resize(3, 3);
    dot_W.resize(3);
    dot_W.clear();
    H.resize(3, 3);
    H.clear();
    invH.resize(3, 3);
    invH.clear();
    dot_u.resize(3);
    for (int dd = 0; dd < 3; dd++) {
      H(dd, dd) = 1;
      invH(dd, dd) = 1;
    }
 
    a_res.resize(3);
    int nb_gauss_pts = row_data.getN().size1();
    commonData.valVel.resize(nb_gauss_pts);
    commonData.jacVelRowPtr.resize(nb_gauss_pts);
    commonData.jacVel.resize(nb_gauss_pts);
 
    int nb_active_vars = 0;
    for (int gg = 0; gg < nb_gauss_pts; gg++) {
 
      if (gg == 0) {
 
        trace_on(tAg);
 
        for (int nn1 = 0; nn1 < 3; nn1++) { // 0
          v[nn1] <<=
              commonData.dataAtGaussPts[commonData.spatialVelocities][gg][nn1];
          nb_active_vars++;
        }
        for (int nn1 = 0; nn1 < 3; nn1++) { // 3
          dot_w[nn1] <<=
              commonData.dataAtGaussPts["DOT_" + commonData.spatialPositions]
                                       [gg][nn1];
          nb_active_vars++;
        }
        if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions]
                .size() > 0) {
          for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3 = 6
            for (int nn2 = 0; nn2 < 3; nn2++) {
              h(nn1, nn2) <<=
                  commonData.gradAtGaussPts[commonData.spatialPositions][gg](
                      nn1, nn2);
              if (fieldDisp) {
                if (nn1 == nn2) {
                  h(nn1, nn2) += 1;
                }
              }
              nb_active_vars++;
            }
          }
          for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
            dot_W[nn1] <<=
                commonData
                    .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
            nb_active_vars++;
          }
        }
        if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
          for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+3
            for (int nn2 = 0; nn2 < 3; nn2++) {
              H(nn1, nn2) <<=
                  commonData.gradAtGaussPts[commonData.meshPositions][gg](nn1,
                                                                          nn2);
              nb_active_vars++;
            }
          }
        }
        detH = 1;
 
        FTensor::Index<'i', 3> i;
        FTensor::Index<'j', 3> j;
        FTensor::Index<'k', 3> k;
 
        auto t_F = getFTensor2FromArray3by3(F, FTensor::Number<0>(), 0);
        auto t_h = getFTensor2FromArray3by3(h, FTensor::Number<0>(), 0);
        auto t_H = getFTensor2FromArray3by3(H, FTensor::Number<0>(), 0);
        auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
        auto t_dot_u =
            FTensor::Tensor1<adouble *, 3>{&dot_u[0], &dot_u[1], &dot_u[2]};
        auto t_dot_w =
            FTensor::Tensor1<adouble *, 3>{&dot_w[0], &dot_w[1], &dot_w[2]};
        auto t_dot_W =
            FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
        auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
        auto t_a_res =
            FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
 
        if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
          detH = determinantTensor3by3(H);
          CHKERR invertTensor3by3(H, detH, invH);
          t_F(i, j) = t_h(i, k) * t_invH(k, j);
        } else {
          t_F(i, j) = t_h(i, j);
        }
 
        t_dot_u(i) = t_dot_w(i) + t_F(i, j) * t_dot_W(j);
        t_a_res(i) = t_v(i) - t_dot_u(i);
        t_a_res(i) *= detH;
 
        // dependant
        VectorDouble &res = commonData.valVel[gg];
        res.resize(3);
        for (int rr = 0; rr < 3; rr++) {
          a_res[rr] >>= res[rr];
        }
        trace_off();
      }
 
      active.resize(nb_active_vars);
      int aa = 0;
      for (int nn1 = 0; nn1 < 3; nn1++) {
        active[aa++] =
            commonData.dataAtGaussPts[commonData.spatialVelocities][gg][nn1];
      }
      for (int nn1 = 0; nn1 < 3; nn1++) {
        active[aa++] =
            commonData
                .dataAtGaussPts["DOT_" + commonData.spatialPositions][gg][nn1];
      }
      if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() >
          0) {
        for (int nn1 = 0; nn1 < 3; nn1++) {
          for (int nn2 = 0; nn2 < 3; nn2++) {
            if (fieldDisp && nn1 == nn2) {
              active[aa++] =
                  commonData.gradAtGaussPts[commonData.spatialPositions][gg](
                      nn1, nn2) +
                  1;
            } else {
              active[aa++] =
                  commonData.gradAtGaussPts[commonData.spatialPositions][gg](
                      nn1, nn2);
            }
          }
        }
        for (int nn1 = 0; nn1 < 3; nn1++) {
          active[aa++] =
              commonData
                  .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
        }
      }
      if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
        for (int nn1 = 0; nn1 < 3; nn1++) {
          for (int nn2 = 0; nn2 < 3; nn2++) {
            active[aa++] =
                commonData.gradAtGaussPts[commonData.meshPositions][gg](nn1,
                                                                        nn2);
          }
        }
      }
 
      if (!jAcobian) {
        VectorDouble &res = commonData.valVel[gg];
        if (gg > 0) {
          res.resize(3);
          int r;
          r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
          if (r != 3) {
            SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                    "ADOL-C function evaluation with error");
          }
        }
        double val = getVolume() * getGaussPts()(3, gg);
        res *= val;
      } else {
        commonData.jacVelRowPtr[gg].resize(3);
        commonData.jacVel[gg].resize(3, nb_active_vars);
        for (int nn1 = 0; nn1 < 3; nn1++) {
          (commonData.jacVelRowPtr[gg])[nn1] = &(commonData.jacVel[gg](nn1, 0));
        }
        int r;
        r = jacobian(tAg, 3, nb_active_vars, &active[0],
                     &(commonData.jacVelRowPtr[gg])[0]);
        if (r != 3) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
                  "ADOL-C function evaluation with error");
        }
        double val = getVolume() * getGaussPts()(3, gg);
        commonData.jacVel[gg] *= val;
        // std::cerr << gg << " : " << commonData.jacVel[gg] << std::endl;
      }
    }
  }
  MoFEMFunctionReturnHot(0);
}

Member Data Documentation

a_res

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

Definition at line 291 of file ConvectiveMassElement.hpp.

active

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

Definition at line 295 of file ConvectiveMassElement.hpp.

commonData

CommonData& ConvectiveMassElement::OpVelocityJacobian::commonData

Definition at line 284 of file ConvectiveMassElement.hpp.

dAta

BlockData& ConvectiveMassElement::OpVelocityJacobian::dAta

Definition at line 283 of file ConvectiveMassElement.hpp.

detH

adouble ConvectiveMassElement::OpVelocityJacobian::detH

Definition at line 293 of file ConvectiveMassElement.hpp.

dot_u

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpVelocityJacobian::dot_u

Definition at line 291 of file ConvectiveMassElement.hpp.

dot_w

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpVelocityJacobian::dot_w

Definition at line 291 of file ConvectiveMassElement.hpp.

dot_W

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

Definition at line 291 of file ConvectiveMassElement.hpp.

F

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

Definition at line 292 of file ConvectiveMassElement.hpp.

fieldDisp

bool ConvectiveMassElement::OpVelocityJacobian::fieldDisp

Definition at line 286 of file ConvectiveMassElement.hpp.

h

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

Definition at line 292 of file ConvectiveMassElement.hpp.

H

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

Definition at line 292 of file ConvectiveMassElement.hpp.

invH

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

Definition at line 292 of file ConvectiveMassElement.hpp.

jAcobian

bool ConvectiveMassElement::OpVelocityJacobian::jAcobian

Definition at line 286 of file ConvectiveMassElement.hpp.

tAg

int ConvectiveMassElement::OpVelocityJacobian::tAg

Definition at line 285 of file ConvectiveMassElement.hpp.

v

VectorBoundedArray<adouble, 3> ConvectiveMassElement::OpVelocityJacobian::v

Definition at line 291 of file ConvectiveMassElement.hpp.

---

The documentation for this struct was generated from the following files:

• users_modules/basic_finite_elements/src/ConvectiveMassElement.hpp
• users_modules/basic_finite_elements/src/impl/ConvectiveMassElement.cpp