v0.15.0
Loading...
Searching...
No Matches
Functions | Variables
save_tensor_data.cpp File Reference

Generate and save on tags thermal stress according to SSLV116 test. More...

#include <MoFEM.hpp>

Go to the source code of this file.

Functions

auto voigt_notation ()
 
auto tensor_notation ()
 
int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Detailed Description

Generate and save on tags thermal stress according to SSLV116 test.

Definition in file save_tensor_data.cpp.

Function Documentation

◆ main()

int main ( int argc,
char * argv[] )

Definition at line 45 of file save_tensor_data.cpp.

45 {
46 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
47
48 try {
49
50 moab::Core mb_instance;
51 moab::Interface &moab = mb_instance;
52
53 // Create MoFEM database
54 MoFEM::Core core(moab);
55 MoFEM::Interface &m_field = core;
56
57 char mesh_file[255] = "mesh.h5m";
58 char output_file[255] = "out.h5m";
59 char tag_name[255] = "INTERNAL_STRESS";
60
61 int tag_len = 9;
62
63 PetscBool use_voigt_notation = PETSC_FALSE;
64
65 PetscOptionsBegin(m_field.get_comm(), "", "save thermal stress test",
66 "none");
67 CHKERR PetscOptionsString("-file_name", " mesh file name", "", "mesh.h5m",
68 mesh_file, 255, PETSC_NULLPTR);
69 CHKERR PetscOptionsString("-output_file", "output file name", "", "out.h5m",
70 output_file, 255, PETSC_NULLPTR);
71 CHKERR PetscOptionsString("-tag_name", "Stress tag name", "",
72 "INTERNAL_STRESS", tag_name, 255, PETSC_NULLPTR);
73 CHKERR PetscOptionsBool("-voigt_notation", "use voigt notation", "",
74 use_voigt_notation, &use_voigt_notation,
75 PETSC_NULLPTR);
76
77 PetscOptionsEnd();
78
79 MOFEM_LOG_CHANNEL("WORLD");
80 MOFEM_LOG_TAG("WORLD", "save_tensor_data");
81
82 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
83 auto moab_comm_wrap =
84 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
85 if (pcomm == NULL)
86 pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
87
88 const std::string options = "PARALLEL=READ_PART;"
89 "PARALLEL_RESOLVE_SHARED_ENTS;"
90 "PARTITION=PARALLEL_PARTITION;";
91 CHKERR moab.load_file(mesh_file, 0, options.c_str());
92
93 Range src_elems;
94 CHKERR pcomm->get_part_entities(src_elems, 3);
95
96 Tag data_tag;
97 double def_val[9];
98 bzero(def_val, 9 * sizeof(double));
99 CHKERR moab.tag_get_handle(tag_name, tag_len, MB_TYPE_DOUBLE, data_tag,
100 MB_TAG_CREAT | MB_TAG_DENSE, &def_val);
101
102 // Set data on the mesh for SSLV116 test
103 std::vector<double> spos;
104 spos.resize(3 * src_elems.size());
105 CHKERR moab.get_coords(src_elems, &spos[0]);
106
107 VectorDouble data(9 * src_elems.size(), 0.0);
108 auto t_data = getFTensor1FromArray<9, 9>(data);
109
110 // temperature distribution for SSLV116 test
111 // see https://code-aster.org/V2/doc/v14/en/man_v/v3/v3.04.116.pdf
112 auto temperature = [](double z) -> double {
113 return (z < -1) ? (120.0 / -9.0) * (z + 10) + 250
114 : (z <= 2) ? (50.0 / -3.0) * (z + 1) + 130
115 : (170.0 / 8.0) * (z - 2) + 80;
116 };
117
118 // material properties for SSLV116 test
119 double young_modulus = 2e11;
120 double poisson_ratio = 0.3;
121 double alpha = 1e-5;
122 double temp_0 = 250;
123
124 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
125 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
126
127 FTensor::Index<'i', 3> i;
128 FTensor::Index<'j', 3> j;
129 FTensor::Index<'k', 3> k;
130 FTensor::Index<'l', 3> l;
131
132 // computing elasticity matrix
133 MatrixDouble mat_D_ptr;
134 mat_D_ptr.resize(6, 6, false);
135 auto t_D = getFTensor4DdgFromPtr<3, 3, 0>(&*(mat_D_ptr.data().begin()));
136 constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
137 t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
138 (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
139 t_kd(i, j) * t_kd(k, l);
140
141 FTensor::Index<'L', 9> L;
142 FTensor::Dg<double, 3, 9> t_L;
143 if (use_voigt_notation) {
144 t_L = voigt_notation();
145 } else {
146 t_L = tensor_notation();
147 }
148
149 FTensor::Tensor2<double, 3, 3> t_stress;
150
151 bool tets_only = true;
152
153 // computing thermal stress field (const on each element)
154 int ielem = 0;
155 for (auto sit = src_elems.begin(); sit != src_elems.end(); sit++) {
156 double z = spos[3 * ielem + 2];
157 double temp = temperature(z);
158 t_stress(i, j) = -t_D(i, j, k, l) * t_kd(k, l) * (temp - temp_0) * alpha;
159 t_data(L) = t_L(i, j, L) * t_stress(i, j);
160 ++t_data;
161 ++ielem;
162 if (tets_only && type_from_handle(*sit) != MBTET) {
163 tets_only = false;
164 }
165 }
166
167 CHKERR moab.tag_set_data(data_tag, src_elems, &data[0]);
168
169 // Check the integrated stress for the SSLV116 test
170 if (tets_only) {
171 std::vector<double> data_integ(tag_len, 0.0);
172 for (auto &tet : src_elems) {
173
174 std::vector<double> tet_data(tag_len, 0.0);
175 CHKERR moab.tag_get_data(data_tag, &tet, 1, &tet_data[0]);
176
177 const EntityHandle *vert_conn;
178 int vert_num;
179 CHKERR moab.get_connectivity(tet, vert_conn, vert_num, true);
180 std::vector<double> vpos(3 * vert_num);
181 CHKERR moab.get_coords(vert_conn, vert_num, &vpos[0]);
182 double vol = Tools::tetVolume(vpos.data());
183
184 for (int itag = 0; itag < tag_len; itag++) {
185 data_integ[itag] += tet_data[itag] * vol;
186 }
187 }
188
189 std::vector<double> global_data_integ(tag_len, 0.0);
190 MPI_Allreduce(&data_integ[0], &global_data_integ[0], tag_len, MPI_DOUBLE,
191 MPI_SUM, m_field.get_comm());
192
193 std::string temp = "Integrated stress for sslv116 test: ";
194 for (int i = 0; i < 9; ++i) {
195 std::stringstream ss;
196 ss << std::scientific << std::setprecision(3) << global_data_integ[i];
197 temp += ss.str() + " ";
198 }
199 MOFEM_LOG("WORLD", Sev::inform) << temp;
200
201 double non_zero_val = 1.655e12;
202 double non_zero_tol = 2e-3;
203 std::vector<int> non_zero_inds(3);
204 std::vector<int> zero_inds(6);
205
206 if (use_voigt_notation) {
207 non_zero_inds = {0, 1, 2};
208 zero_inds = {3, 4, 5, 6, 7, 8};
209 } else {
210 non_zero_inds = {0, 4, 8};
211 zero_inds = {1, 2, 3, 5, 6, 7};
212 }
213 bool non_zero_check =
214 all_of(begin(non_zero_inds), end(non_zero_inds), [&](int i) {
215 return abs(global_data_integ[i] - non_zero_val) / non_zero_val <
216 non_zero_tol;
217 });
218 bool zero_check = all_of(begin(zero_inds), end(zero_inds), [&](int i) {
219 return abs(global_data_integ[i]) < 1e-12;
220 });
221 if (!non_zero_check || !zero_check) {
222 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
223 "Wrong value of the integrated stress");
224 }
225 } else {
226 MOFEM_LOG("WORLD", Sev::warning)
227 << "Mesh has non-tetrahedral elements, skipping the test based on "
228 "integrating the stress in the volume";
229 }
230
231 CHKERR moab.write_file(output_file, "MOAB", "PARALLEL=WRITE_PART");
232 }
233 CATCH_ERRORS;
234
235 MoFEM::Core::Finalize();
236
237 return 0;
238}
FTensor::Kronecker_Delta_symmetric
Kronecker Delta class symmetric.
Definition Kronecker_Delta.hpp:49
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
bulk_modulus_K
double bulk_modulus_K
Definition dynamic_first_order_con_law.cpp:95
shear_modulus_G
double shear_modulus_G
Definition dynamic_first_order_con_law.cpp:96
t_kd
constexpr auto t_kd
Definition free_surface.cpp:139
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition LogManager.hpp:347
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition LogManager.hpp:292
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
l
FTensor::Index< 'l', 3 > l
Definition matrix_function.cpp:21
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition Templates.hpp:1932
MoFEM::getFTensor4DdgFromPtr
static auto getFTensor4DdgFromPtr(T *ptr)
Definition Templates.hpp:447
MoFEM::getFTensor1FromArray
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
Definition Templates.hpp:1209
young_modulus
double young_modulus
Young modulus.
Definition plastic.cpp:125
poisson_ratio
double poisson_ratio
Poisson ratio.
Definition plastic.cpp:126
help
static char help[]
Definition save_tensor_data.cpp:43
tensor_notation
auto tensor_notation()
Definition save_tensor_data.cpp:25
voigt_notation
auto voigt_notation()
Definition save_tensor_data.cpp:10
temp
void temp(int x, int y=10)
Definition simple.cpp:4
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:118
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::Tools::tetVolume
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition Tools.cpp:30

◆ tensor_notation()

auto tensor_notation ( )
inline

Definition at line 25 of file save_tensor_data.cpp.

25 {
26 FTensor::Index<'i', 3> i;
27 FTensor::Index<'j', 3> j;
28 FTensor::Index<'L', 9> L;
29 FTensor::Dg<double, 3, 9> t_L;
30 t_L(i, j, L) = 0;
31 t_L(0, 0, 0) = 1;
32 t_L(0, 1, 1) = 1;
33 t_L(0, 2, 2) = 1;
34 t_L(1, 0, 3) = 1;
35 t_L(1, 1, 4) = 1;
36 t_L(1, 2, 5) = 1;
37 t_L(2, 0, 6) = 1;
38 t_L(2, 1, 7) = 1;
39 t_L(2, 2, 8) = 1;
40 return t_L;
41}

◆ voigt_notation()

auto voigt_notation ( )
inline

Definition at line 10 of file save_tensor_data.cpp.

10 {
11 FTensor::Index<'i', 3> i;
12 FTensor::Index<'j', 3> j;
13 FTensor::Index<'L', 9> L;
14 FTensor::Dg<double, 3, 9> t_L;
15 t_L(i, j, L) = 0;
16 t_L(0, 0, 0) = 1;
17 t_L(1, 1, 1) = 1;
18 t_L(2, 2, 2) = 1;
19 t_L(0, 1, 3) = 1;
20 t_L(0, 2, 4) = 1;
21 t_L(1, 2, 5) = 1;
22 return t_L;
23}

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 43 of file save_tensor_data.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow