v0.16.0
Loading...
Searching...
No Matches
Public Member Functions | Public Attributes | List of all members
ConvectiveMassElement::OpVelocityJacobian Struct Reference
User modules » Mass Element

#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

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

793 : VolumeElementForcesAndSourcesCore::UserDataOperator(
794 field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
795 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
796 fieldDisp(false) {}
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROW
@ OPROW
operator doWork function is executed on FE rows
Definition ForcesAndSourcesCore.hpp:602

Member Function Documentation

◆ doWork()

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

Definition at line 798 of file ConvectiveMassElement.cpp.

799 {
800 MoFEMFunctionBeginHot;
801
802 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
803 dAta.tEts.end()) {
804 MoFEMFunctionReturnHot(0);
805 }
806
807 // do it only once, no need to repeat this for edges,faces or tets
808 if (row_type != MBVERTEX)
809 MoFEMFunctionReturnHot(0);
810
811 int nb_dofs = row_data.getIndices().size();
812 if (nb_dofs == 0)
813 MoFEMFunctionReturnHot(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] <<=
849 commonData.dataAtGaussPts[commonData.spatialVelocities][gg][nn1];
850 nb_active_vars++;
851 }
852 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
853 dot_w[nn1] <<=
854 commonData.dataAtGaussPts["DOT_" + commonData.spatialPositions]
855 [gg][nn1];
856 nb_active_vars++;
857 }
858 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions]
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) <<=
863 commonData.gradAtGaussPts[commonData.spatialPositions][gg](
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] <<=
875 commonData
876 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
877 nb_active_vars++;
878 }
879 }
880 if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
881 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+3
882 for (int nn2 = 0; nn2 < 3; nn2++) {
883 H(nn1, nn2) <<=
884 commonData.gradAtGaussPts[commonData.meshPositions][gg](nn1,
885 nn2);
886 nb_active_vars++;
887 }
888 }
889 }
890 detH = 1;
891
892 FTensor::Index<'i', 3> i;
893 FTensor::Index<'j', 3> j;
894 FTensor::Index<'k', 3> k;
895
896 using ActiveMat3 = MatrixBoundedArray<adouble, 9>;
897 auto t_F = GetFTensor2FromArrayImpl<3, 3, 0, ActiveMat3>::get(F, 0, 0);
898 auto t_h = GetFTensor2FromArrayImpl<3, 3, 0, ActiveMat3>::get(h, 0, 0);
899 auto t_H = GetFTensor2FromArrayImpl<3, 3, 0, ActiveMat3>::get(H, 0, 0);
900 auto t_invH =
901 GetFTensor2FromArrayImpl<3, 3, 0, ActiveMat3>::get(invH, 0, 0);
902 auto t_dot_u =
903 FTensor::Tensor1<adouble *, 3>{&dot_u[0], &dot_u[1], &dot_u[2]};
904 auto t_dot_w =
905 FTensor::Tensor1<adouble *, 3>{&dot_w[0], &dot_w[1], &dot_w[2]};
906 auto t_dot_W =
907 FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
908 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
909 auto t_a_res =
910 FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
911
912 if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
913 detH = determinantTensor3by3(H);
914 CHKERR invertTensor3by3(H, detH, invH);
915 t_F(i, j) = t_h(i, k) * t_invH(k, j);
916 } else {
917 t_F(i, j) = t_h(i, j);
918 }
919
920 t_dot_u(i) = t_dot_w(i) + t_F(i, j) * t_dot_W(j);
921 t_a_res(i) = t_v(i) - t_dot_u(i);
922 t_a_res(i) *= detH;
923
924 // dependant
925 VectorDouble &res = commonData.valVel[gg];
926 res.resize(3);
927 for (int rr = 0; rr < 3; rr++) {
928 a_res[rr] >>= res[rr];
929 }
930 trace_off();
931 }
932
933 active.resize(nb_active_vars);
934 int aa = 0;
935 for (int nn1 = 0; nn1 < 3; nn1++) {
936 active[aa++] =
937 commonData.dataAtGaussPts[commonData.spatialVelocities][gg][nn1];
938 }
939 for (int nn1 = 0; nn1 < 3; nn1++) {
940 active[aa++] =
941 commonData
942 .dataAtGaussPts["DOT_" + commonData.spatialPositions][gg][nn1];
943 }
944 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() >
945 0) {
946 for (int nn1 = 0; nn1 < 3; nn1++) {
947 for (int nn2 = 0; nn2 < 3; nn2++) {
948 if (fieldDisp && nn1 == nn2) {
949 active[aa++] =
950 commonData.gradAtGaussPts[commonData.spatialPositions][gg](
951 nn1, nn2) +
952 1;
953 } else {
954 active[aa++] =
955 commonData.gradAtGaussPts[commonData.spatialPositions][gg](
956 nn1, nn2);
957 }
958 }
959 }
960 for (int nn1 = 0; nn1 < 3; nn1++) {
961 active[aa++] =
962 commonData
963 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
964 }
965 }
966 if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
967 for (int nn1 = 0; nn1 < 3; nn1++) {
968 for (int nn2 = 0; nn2 < 3; nn2++) {
969 active[aa++] =
970 commonData.gradAtGaussPts[commonData.meshPositions][gg](nn1,
971 nn2);
972 }
973 }
974 }
975
976 if (!jAcobian) {
977 VectorDouble &res = commonData.valVel[gg];
978 if (gg > 0) {
979 res.resize(3);
980 int r;
981 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
982 if (r != 3) {
983 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
984 "ADOL-C function evaluation with error");
985 }
986 }
987 double val = getVolume() * getGaussPts()(3, gg);
988 res *= val;
989 } else {
990 commonData.jacVelRowPtr[gg].resize(3);
991 commonData.jacVel[gg].resize(3, nb_active_vars);
992 for (int nn1 = 0; nn1 < 3; nn1++) {
993 (commonData.jacVelRowPtr[gg])[nn1] = &(commonData.jacVel[gg](nn1, 0));
994 }
995 int r;
996 r = jacobian(tAg, 3, nb_active_vars, &active[0],
997 &(commonData.jacVelRowPtr[gg])[0]);
998 if (r != 3) {
999 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
1000 "ADOL-C function evaluation with error");
1001 }
1002 double val = getVolume() * getGaussPts()(3, gg);
1003 commonData.jacVel[gg] *= val;
1004 // std::cerr << gg << " : " << commonData.jacVel[gg] << std::endl;
1005 }
1006 }
1007 }
1008 MoFEMFunctionReturnHot(0);
1009}
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MOFEM_OPERATION_UNSUCCESSFUL
@ MOFEM_OPERATION_UNSUCCESSFUL
Definition definitions.h:34
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:22
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::MatrixBoundedArray
ublas::matrix< T, ublas::row_major, ublas::bounded_array< T, N > > MatrixBoundedArray
Definition Types.hpp:104
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:2276
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition Templates.hpp:2228
sdf.r
int r
Definition sdf.py:205
ConvectiveMassElement::CommonData::dataAtGaussPts
std::map< std::string, std::vector< VectorDouble > > dataAtGaussPts
Definition ConvectiveMassElement.hpp:133
ConvectiveMassElement::CommonData::jacVelRowPtr
std::vector< std::vector< double * > > jacVelRowPtr
Definition ConvectiveMassElement.hpp:139
ConvectiveMassElement::CommonData::gradAtGaussPts
std::map< std::string, std::vector< MatrixDouble > > gradAtGaussPts
Definition ConvectiveMassElement.hpp:134
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:1600
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of degrees of freedom on entity.
Definition EntitiesFieldData.hpp:1523
MoFEM::GetFTensor2FromArrayImpl
Get FTensor2 from array.
Definition Templates.hpp:660

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:
Generated by Doxygen 1.9.8 and hosted at University of Glasgow