v0.13.1
Public Member Functions | Public Attributes | List of all members
ConvectiveMassElement::OpVelocityJacobian Struct Reference

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

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

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

BlockDatadAta
 
CommonDatacommonData
 
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< doubleactive
 

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

Member Function Documentation

◆ doWork()

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

Definition at line 797 of file ConvectiveMassElement.cpp.

799 {
801
802 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
803 dAta.tEts.end()) {
805 }
806
807 // do it only once, no need to repeat this for edges,faces or tets
808 if (row_type != MBVERTEX)
810
811 int nb_dofs = row_data.getIndices().size();
812 if (nb_dofs == 0)
814
815 {
816
817 v.resize(3);
818 dot_w.resize(3);
819 h.resize(3, 3);
820 h.clear();
821 F.resize(3, 3);
822 dot_W.resize(3);
823 dot_W.clear();
824 H.resize(3, 3);
825 H.clear();
826 invH.resize(3, 3);
827 invH.clear();
828 dot_u.resize(3);
829 for (int dd = 0; dd < 3; dd++) {
830 H(dd, dd) = 1;
831 invH(dd, dd) = 1;
832 }
833
834 a_res.resize(3);
835 int nb_gauss_pts = row_data.getN().size1();
836 commonData.valVel.resize(nb_gauss_pts);
837 commonData.jacVelRowPtr.resize(nb_gauss_pts);
838 commonData.jacVel.resize(nb_gauss_pts);
839
840 int nb_active_vars = 0;
841 for (int gg = 0; gg < nb_gauss_pts; gg++) {
842
843 if (gg == 0) {
844
845 trace_on(tAg);
846
847 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
848 v[nn1] <<=
850 nb_active_vars++;
851 }
852 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
853 dot_w[nn1] <<=
855 [gg][nn1];
856 nb_active_vars++;
857 }
859 .size() > 0) {
860 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3 = 6
861 for (int nn2 = 0; nn2 < 3; nn2++) {
862 h(nn1, nn2) <<=
864 nn1, nn2);
865 if (fieldDisp) {
866 if (nn1 == nn2) {
867 h(nn1, nn2) += 1;
868 }
869 }
870 nb_active_vars++;
871 }
872 }
873 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
874 dot_W[nn1] <<=
876 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
877 nb_active_vars++;
878 }
879 }
881 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+3
882 for (int nn2 = 0; nn2 < 3; nn2++) {
883 H(nn1, nn2) <<=
885 nn2);
886 nb_active_vars++;
887 }
888 }
889 }
890 detH = 1;
891
895
900 auto t_dot_u =
902 auto t_dot_w =
904 auto t_dot_W =
906 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
907 auto t_a_res =
909
913 t_F(i, j) = t_h(i, k) * t_invH(k, j);
914 } else {
915 t_F(i, j) = t_h(i, j);
916 }
917
918 t_dot_u(i) = t_dot_w(i) + t_F(i, j) * t_dot_W(j);
919 t_a_res(i) = t_v(i) - t_dot_u(i);
920 t_a_res(i) *= detH;
921
922 // dependant
923 VectorDouble &res = commonData.valVel[gg];
924 res.resize(3);
925 for (int rr = 0; rr < 3; rr++) {
926 a_res[rr] >>= res[rr];
927 }
928 trace_off();
929 }
930
931 active.resize(nb_active_vars);
932 int aa = 0;
933 for (int nn1 = 0; nn1 < 3; nn1++) {
934 active[aa++] =
936 }
937 for (int nn1 = 0; nn1 < 3; nn1++) {
938 active[aa++] =
940 .dataAtGaussPts["DOT_" + commonData.spatialPositions][gg][nn1];
941 }
943 0) {
944 for (int nn1 = 0; nn1 < 3; nn1++) {
945 for (int nn2 = 0; nn2 < 3; nn2++) {
946 if (fieldDisp && nn1 == nn2) {
947 active[aa++] =
949 nn1, nn2) +
950 1;
951 } else {
952 active[aa++] =
954 nn1, nn2);
955 }
956 }
957 }
958 for (int nn1 = 0; nn1 < 3; nn1++) {
959 active[aa++] =
961 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
962 }
963 }
965 for (int nn1 = 0; nn1 < 3; nn1++) {
966 for (int nn2 = 0; nn2 < 3; nn2++) {
967 active[aa++] =
969 nn2);
970 }
971 }
972 }
973
974 if (!jAcobian) {
975 VectorDouble &res = commonData.valVel[gg];
976 if (gg > 0) {
977 res.resize(3);
978 int r;
979 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
980 if (r != 3) {
981 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
982 "ADOL-C function evaluation with error");
983 }
984 }
985 double val = getVolume() * getGaussPts()(3, gg);
986 res *= val;
987 } else {
988 commonData.jacVelRowPtr[gg].resize(3);
989 commonData.jacVel[gg].resize(3, nb_active_vars);
990 for (int nn1 = 0; nn1 < 3; nn1++) {
991 (commonData.jacVelRowPtr[gg])[nn1] = &(commonData.jacVel[gg](nn1, 0));
992 }
993 int r;
994 r = jacobian(tAg, 3, nb_active_vars, &active[0],
995 &(commonData.jacVelRowPtr[gg])[0]);
996 if (r != 3) {
997 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
998 "ADOL-C function evaluation with error");
999 }
1000 double val = getVolume() * getGaussPts()(3, gg);
1001 commonData.jacVel[gg] *= val;
1002 // std::cerr << gg << " : " << commonData.jacVel[gg] << std::endl;
1003 }
1004 }
1005 }
1007}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition: definitions.h:34
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
auto getFTensor2FromArray3by3(ublas::matrix< T, L, A > &data, const FTensor::Number< S > &, const size_t rr)
Definition: Templates.hpp:972
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
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1193
const double r
rate factor
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
std::vector< std::vector< double * > > jacVelRowPtr
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.

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: