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

792 : VolumeElementForcesAndSourcesCore::UserDataOperator(
793 field_name, ForcesAndSourcesCore::UserDataOperator::OPROW),
794 dAta(data), commonData(common_data), tAg(tag), jAcobian(jacobian),
795 fieldDisp(false) {}
field_name
constexpr auto field_name
Definition poisson_2d_homogeneous.cpp:13

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.

798 {
799 MoFEMFunctionBeginHot;
800
801 if (dAta.tEts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
802 dAta.tEts.end()) {
803 MoFEMFunctionReturnHot(0);
804 }
805
806 // do it only once, no need to repeat this for edges,faces or tets
807 if (row_type != MBVERTEX)
808 MoFEMFunctionReturnHot(0);
809
810 int nb_dofs = row_data.getIndices().size();
811 if (nb_dofs == 0)
812 MoFEMFunctionReturnHot(0);
813
814 {
815
816 v.resize(3);
817 dot_w.resize(3);
818 h.resize(3, 3);
819 h.clear();
820 F.resize(3, 3);
821 dot_W.resize(3);
822 dot_W.clear();
823 H.resize(3, 3);
824 H.clear();
825 invH.resize(3, 3);
826 invH.clear();
827 dot_u.resize(3);
828 for (int dd = 0; dd < 3; dd++) {
829 H(dd, dd) = 1;
830 invH(dd, dd) = 1;
831 }
832
833 a_res.resize(3);
834 int nb_gauss_pts = row_data.getN().size1();
835 commonData.valVel.resize(nb_gauss_pts);
836 commonData.jacVelRowPtr.resize(nb_gauss_pts);
837 commonData.jacVel.resize(nb_gauss_pts);
838
839 int nb_active_vars = 0;
840 for (int gg = 0; gg < nb_gauss_pts; gg++) {
841
842 if (gg == 0) {
843
844 trace_on(tAg);
845
846 for (int nn1 = 0; nn1 < 3; nn1++) { // 0
847 v[nn1] <<=
848 commonData.dataAtGaussPts[commonData.spatialVelocities][gg][nn1];
849 nb_active_vars++;
850 }
851 for (int nn1 = 0; nn1 < 3; nn1++) { // 3
852 dot_w[nn1] <<=
853 commonData.dataAtGaussPts["DOT_" + commonData.spatialPositions]
854 [gg][nn1];
855 nb_active_vars++;
856 }
857 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions]
858 .size() > 0) {
859 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3 = 6
860 for (int nn2 = 0; nn2 < 3; nn2++) {
861 h(nn1, nn2) <<=
862 commonData.gradAtGaussPts[commonData.spatialPositions][gg](
863 nn1, nn2);
864 if (fieldDisp) {
865 if (nn1 == nn2) {
866 h(nn1, nn2) += 1;
867 }
868 }
869 nb_active_vars++;
870 }
871 }
872 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9
873 dot_W[nn1] <<=
874 commonData
875 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
876 nb_active_vars++;
877 }
878 }
879 if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
880 for (int nn1 = 0; nn1 < 3; nn1++) { // 3+3+9+3
881 for (int nn2 = 0; nn2 < 3; nn2++) {
882 H(nn1, nn2) <<=
883 commonData.gradAtGaussPts[commonData.meshPositions][gg](nn1,
884 nn2);
885 nb_active_vars++;
886 }
887 }
888 }
889 detH = 1;
890
891 FTensor::Index<'i', 3> i;
892 FTensor::Index<'j', 3> j;
893 FTensor::Index<'k', 3> k;
894
895 auto t_F = getFTensor2FromArray3by3(F, FTensor::Number<0>(), 0);
896 auto t_h = getFTensor2FromArray3by3(h, FTensor::Number<0>(), 0);
897 auto t_H = getFTensor2FromArray3by3(H, FTensor::Number<0>(), 0);
898 auto t_invH = getFTensor2FromArray3by3(invH, FTensor::Number<0>(), 0);
899 auto t_dot_u =
900 FTensor::Tensor1<adouble *, 3>{&dot_u[0], &dot_u[1], &dot_u[2]};
901 auto t_dot_w =
902 FTensor::Tensor1<adouble *, 3>{&dot_w[0], &dot_w[1], &dot_w[2]};
903 auto t_dot_W =
904 FTensor::Tensor1<adouble *, 3>{&dot_W[0], &dot_W[1], &dot_W[2]};
905 auto t_v = FTensor::Tensor1<adouble *, 3>{&v[0], &v[1], &v[2]};
906 auto t_a_res =
907 FTensor::Tensor1<adouble *, 3>{&a_res[0], &a_res[1], &a_res[2]};
908
909 if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
910 detH = determinantTensor3by3(H);
911 CHKERR invertTensor3by3(H, detH, invH);
912 t_F(i, j) = t_h(i, k) * t_invH(k, j);
913 } else {
914 t_F(i, j) = t_h(i, j);
915 }
916
917 t_dot_u(i) = t_dot_w(i) + t_F(i, j) * t_dot_W(j);
918 t_a_res(i) = t_v(i) - t_dot_u(i);
919 t_a_res(i) *= detH;
920
921 // dependant
922 VectorDouble &res = commonData.valVel[gg];
923 res.resize(3);
924 for (int rr = 0; rr < 3; rr++) {
925 a_res[rr] >>= res[rr];
926 }
927 trace_off();
928 }
929
930 active.resize(nb_active_vars);
931 int aa = 0;
932 for (int nn1 = 0; nn1 < 3; nn1++) {
933 active[aa++] =
934 commonData.dataAtGaussPts[commonData.spatialVelocities][gg][nn1];
935 }
936 for (int nn1 = 0; nn1 < 3; nn1++) {
937 active[aa++] =
938 commonData
939 .dataAtGaussPts["DOT_" + commonData.spatialPositions][gg][nn1];
940 }
941 if (commonData.dataAtGaussPts["DOT_" + commonData.meshPositions].size() >
942 0) {
943 for (int nn1 = 0; nn1 < 3; nn1++) {
944 for (int nn2 = 0; nn2 < 3; nn2++) {
945 if (fieldDisp && nn1 == nn2) {
946 active[aa++] =
947 commonData.gradAtGaussPts[commonData.spatialPositions][gg](
948 nn1, nn2) +
949 1;
950 } else {
951 active[aa++] =
952 commonData.gradAtGaussPts[commonData.spatialPositions][gg](
953 nn1, nn2);
954 }
955 }
956 }
957 for (int nn1 = 0; nn1 < 3; nn1++) {
958 active[aa++] =
959 commonData
960 .dataAtGaussPts["DOT_" + commonData.meshPositions][gg][nn1];
961 }
962 }
963 if (commonData.gradAtGaussPts[commonData.meshPositions].size() > 0) {
964 for (int nn1 = 0; nn1 < 3; nn1++) {
965 for (int nn2 = 0; nn2 < 3; nn2++) {
966 active[aa++] =
967 commonData.gradAtGaussPts[commonData.meshPositions][gg](nn1,
968 nn2);
969 }
970 }
971 }
972
973 if (!jAcobian) {
974 VectorDouble &res = commonData.valVel[gg];
975 if (gg > 0) {
976 res.resize(3);
977 int r;
978 r = ::function(tAg, 3, nb_active_vars, &active[0], &res[0]);
979 if (r != 3) {
980 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
981 "ADOL-C function evaluation with error");
982 }
983 }
984 double val = getVolume() * getGaussPts()(3, gg);
985 res *= val;
986 } else {
987 commonData.jacVelRowPtr[gg].resize(3);
988 commonData.jacVel[gg].resize(3, nb_active_vars);
989 for (int nn1 = 0; nn1 < 3; nn1++) {
990 (commonData.jacVelRowPtr[gg])[nn1] = &(commonData.jacVel[gg](nn1, 0));
991 }
992 int r;
993 r = jacobian(tAg, 3, nb_active_vars, &active[0],
994 &(commonData.jacVelRowPtr[gg])[0]);
995 if (r != 3) {
996 SETERRQ(PETSC_COMM_SELF, MOFEM_OPERATION_UNSUCCESSFUL,
997 "ADOL-C function evaluation with error");
998 }
999 double val = getVolume() * getGaussPts()(3, gg);
1000 commonData.jacVel[gg] *= val;
1001 // std::cerr << gg << " : " << commonData.jacVel[gg] << std::endl;
1002 }
1003 }
1004 }
1005 MoFEMFunctionReturnHot(0);
1006}
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:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
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::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:1622
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:1376
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition Templates.hpp:1574
sdf.r
int r
Definition sdf.py:8
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:1318
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition EntitiesFieldData.hpp:1214

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.12.0 and hosted at University of Glasgow