v0.15.5
Loading...
Searching...
No Matches
Functions | Variables
project_roughness.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Functions

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

Variables

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

Function Documentation

◆ main()

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

Definition at line 28 of file project_roughness.cpp.

28 {
29
30 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
31
32 try {
33
34 moab::Core mb_instance;
35 moab::Interface &moab = mb_instance;
36 int rank;
37 MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
38
39 // Read parameters from line command
40 PetscBool flg_mesh_file = PETSC_TRUE;
41 PetscBool flg_out_file = PETSC_TRUE;
42 PetscBool flg_data_file = PETSC_TRUE;
43
44 char mesh_file_name[255];
45 char out_file_name[255] = "out.h5m";
46 char data_file_name[255];
47
48 int nb_side_nodes = 1;
49 double side_length = 1.0;
50 double mesh_height = 1.0;
51 double mesh_scale = 1.0;
52 double tol = 1e-5;
53 int atom_test = 0;
54
55 PetscOptionsBegin(PETSC_COMM_WORLD, "", "Project Roughness", nullptr);
56 CHKERR PetscOptionsString("-file_name", "mesh file name", "", "mesh.h5m",
57 mesh_file_name, 255, &flg_mesh_file);
58 CHKERR PetscOptionsString("-output_file", "out file name", "", "out.h5m",
59 out_file_name, 255, &flg_out_file);
60 CHKERR PetscOptionsString("-surface_data_file", "data file name", "",
61 "data.txt", data_file_name, 255, &flg_data_file);
62
63 CHKERR PetscOptionsInt("-mesh_side_nodes_num",
64 "number of nodes on each side", "", nb_side_nodes,
65 &nb_side_nodes, PETSC_NULLPTR);
66 CHKERR PetscOptionsReal("-mesh_side_length", "side length", "", side_length,
67 &side_length, PETSC_NULLPTR);
68 CHKERR PetscOptionsReal("-mesh_scale", "mesh scale", "", mesh_scale,
69 &mesh_scale, PETSC_NULLPTR);
70 CHKERR PetscOptionsReal("-mesh_height", "vertical height of the mesh", "",
71 mesh_height, &mesh_height, PETSC_NULLPTR);
72 CHKERR PetscOptionsReal("-tolerance", "Tolerance for coordinate matching",
73 "", tol, &tol, PETSC_NULLPTR);
74 CHKERR PetscOptionsInt("-atom_test", "atom test number", "", 0, &atom_test,
75 PETSC_NULLPTR);
76
77 PetscOptionsEnd();
78
79 MOFEM_LOG_CHANNEL("SELF");
80
81 MOFEM_LOG_C("SELF", Sev::inform, "Mesh side length: %.2f", side_length);
82 MOFEM_LOG_C("SELF", Sev::inform, "Mesh height: %.2f", mesh_height);
83 MOFEM_LOG_C("SELF", Sev::inform, "Number of nodes on each side: %d",
84 nb_side_nodes);
85
86 std::ifstream infile(data_file_name, std::ios::in);
87 std::vector<double> surface_data;
88 std::vector<bool> surface_check;
89
90 if (!infile.is_open()) {
91 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
92 "Could not open surface data file");
93 }
94
95 std::string line;
96 while (std::getline(infile, line)) {
97 std::istringstream iss(line);
98 double num = 0.0;
99 if (iss >> num)
100 surface_data.push_back(num);
101 else
102 continue;
103 }
104 int expected_surface_nodes = nb_side_nodes * nb_side_nodes;
105
106 if ((int)surface_data.size() != expected_surface_nodes) {
107 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
108 "Surface data size mismatch: expected %d nodes (%d x %d), got %d",
109 expected_surface_nodes, nb_side_nodes, nb_side_nodes,
110 (int)surface_data.size());
111 }
112 double max_surf_height =
113 *std::max_element(std::begin(surface_data), std::end(surface_data));
114 for (auto &node_data : surface_data) {
115 node_data -= max_surf_height;
116 }
117
118 surface_check.resize(surface_data.size());
119 std::fill(surface_check.begin(), surface_check.end(), false);
120
121 if (flg_mesh_file != PETSC_TRUE) {
122 SETERRQ(PETSC_COMM_SELF, 1, "*** ERROR -my_file (MESH FILE NEEDED)");
123 }
124
125 // Read mesh to MOAB
126 const char *option;
127 option = ""; //"PARALLEL=BCAST;";//;DEBUG_IO";
128 CHKERR moab.load_file(mesh_file_name, 0, option);
129 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
130 if (pcomm == NULL)
131 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
132
133 MoFEM::Core core(moab);
134 MOFEM_LOG_C("SELF", Sev::inform, "%s",
135 "Projecting surface data on the mesh");
136
137 Range elastic_elems, elastic_nodes;
138 CHKERR moab.get_entities_by_dimension(
139 0, 3, elastic_elems); // Todo: currently for 3D meshes only
140 CHKERR moab.get_connectivity(elastic_elems, elastic_nodes);
141
142 const double edge_length = side_length / (double)(nb_side_nodes - 1);
143
144 MOFEM_LOG_C("SELF", Sev::inform, "Actual edge length per element: %.8f",
145 edge_length);
146
147 auto flat_ind = [&](int i, int j) { return j * nb_side_nodes + i; };
148
149 for (Range::iterator nit = elastic_nodes.begin();
150 nit != elastic_nodes.end(); nit++) {
151 double coords[3];
152 CHKERR moab.get_coords(&*nit, 1, coords);
153 double x = coords[0];
154 double y = coords[1];
155 double z = coords[2];
156 double height_coef = (mesh_height + z) / mesh_height;
157
158 double ksi, eta;
159 double pos_x = x / edge_length;
160 double pos_y = y / edge_length;
161
162 double pos_x_round = std::round(pos_x);
163 double pos_y_round = std::round(pos_y);
164
165 if (std::abs(pos_x_round - pos_x) < tol) {
166 pos_x = pos_x_round;
167 }
168 if (std::abs(pos_y_round - pos_y) < tol) {
169 pos_y = pos_y_round;
170 }
171
172 int i, j, i_next, j_next;
173 i = i_next = (int)std::floor(pos_x);
174 j = j_next = (int)std::floor(pos_y);
175
176 if (std::abs(pos_x - std::floor(pos_x)) < tol) {
177 ksi = -1.0;
178 } else {
179 ksi = (pos_x - std::floor(pos_x)) * 2.0 - 1.0;
180 if (ksi < -1. - tol || ksi > 1. + tol) {
181 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
182 "Wrong ksi for i = %d, j = %d", i, j);
183 }
184 i_next++;
185 }
186 if (std::abs(pos_y - std::floor(pos_y)) < tol) {
187 eta = -1.0;
188 } else {
189 eta = (pos_y - std::floor(pos_y)) * 2.0 - 1.0;
190 if (eta < -1. - tol || eta > 1. + tol) {
191 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
192 "Wrong eta for i = %d, j = %d", i, j);
193 }
194 j_next++;
195 }
196 if (std::abs(z) < tol) {
197 if (std::abs(pos_x - std::floor(pos_x)) > tol ||
198 std::abs(pos_y - std::floor(pos_y)) > tol) {
199 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
200 "Cannot find data for surface node i = %d, j = %d", i, j);
201 } else {
202 if (surface_check[flat_ind(i, j)]) {
203 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
204 "Data used twice for surface node i = %d, j = %d", i, j);
205 } else {
206 surface_check[flat_ind(i, j)] = true;
207 }
208 }
209 }
210 std::array<double, 4> shape_quad;
211 shape_quad[0] = 0.25 * (1.0 - ksi) * (1.0 - eta);
212 shape_quad[1] = 0.25 * (1.0 + ksi) * (1.0 - eta);
213 shape_quad[2] = 0.25 * (1.0 + ksi) * (1.0 + eta);
214 shape_quad[3] = 0.25 * (1.0 - ksi) * (1.0 + eta);
215
216 std::array<double, 4> gap_quad;
217 gap_quad[0] = surface_data[flat_ind(i, j)];
218 gap_quad[1] = surface_data[flat_ind(i_next, j)];
219 gap_quad[2] = surface_data[flat_ind(i_next, j_next)];
220 gap_quad[3] = surface_data[flat_ind(i, j_next)];
221
222 for (int k : {0, 1, 2, 3}) {
223 coords[2] += shape_quad[k] * gap_quad[k] * height_coef;
224 }
225 CHKERR moab.set_coords(&*nit, 1, coords);
226 }
227
228 for (int i = 0; i < nb_side_nodes; i++) {
229 for (int j = 0; j < nb_side_nodes; j++) {
230 if (!surface_check[flat_ind(i, j)]) {
231 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
232 "Data not used for surface node i = %d, j = %d", i, j);
233 }
234 }
235 }
236 MOFEM_LOG_C("SELF", Sev::inform, "%s", "Scaling the mesh");
237
238 for (Range::iterator nit = elastic_nodes.begin();
239 nit != elastic_nodes.end(); nit++) {
240 double coords[3];
241 CHKERR moab.get_coords(&*nit, 1, coords);
242 for (int i : {0, 1, 2}) {
243 coords[i] *= mesh_scale;
244 }
245 CHKERR moab.set_coords(&*nit, 1, coords);
246 }
247
248 Range tets;
249 CHKERR moab.get_entities_by_type(0, MBTET, tets);
250
251 if (tets.empty()) {
252 MOFEM_LOG_C("SELF", Sev::warning, "%s",
253 "No tetrahedral elements found, skipping quality and volume checks");
254 } else {
255 // variables for quality check and volume computation
256 double total_volume = 0.0;
257 double min_quality = std::numeric_limits<double>::max();
258
259 for (auto tet : tets) {
260 const EntityHandle *conn = nullptr;
261 int num_nodes = 0;
262 double coords[12];
263
264 CHKERR moab.get_connectivity(tet, conn, num_nodes, true);
265 CHKERR moab.get_coords(conn, num_nodes, coords);
266
267 double q = Tools::volumeLengthQuality(coords);
268 if (!std::isnormal(q) || q <= 0) {
269 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
270 "invalid tet quality %.6e", q);
271 }
272 min_quality = std::min(min_quality, q);
273
274 double vol = Tools::tetVolume(coords);
275 total_volume += vol;
276 }
277 MOFEM_LOG_C("SELF", Sev::inform, "Min quality of the final mesh = %.6e",
278 min_quality);
279 MOFEM_LOG_C("SELF", Sev::inform, "Total volume of the final mesh: %.4f",
280 total_volume);
281
282 if (atom_test == 1) {
283 double expected_volume = 1.98;
284 double diff = std::abs(total_volume - expected_volume);
285
286 if (diff > tol) {
287 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
288 "Atom test %d failed: Unexpected total tet volume computed: "
289 "%.6e (Expected %.6e)",
290 atom_test, total_volume, expected_volume);
291 }
292 }
293 }
294 MOFEM_LOG_C("SELF", Sev::inform, "Writing output file: %s", out_file_name);
295 EntityHandle rootset = moab.get_root_set();
296 CHKERR moab.write_file(out_file_name, "MOAB", "", &rootset, 1);
297 }
298 CATCH_ERRORS;
299
300 MoFEM::Core::Finalize();
301 return 0;
302}
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition LogManager.hpp:319
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_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition definitions.h:33
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
eta
double eta
Definition free_surface.cpp:192
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
out_file_name
char out_file_name[255]
Definition initial_diffusion.cpp:53
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20
help
static char help[]
Definition project_roughness.cpp:26
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::Tools::tetVolume
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition Tools.cpp:30
MoFEM::Tools::volumeLengthQuality
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition Tools.cpp:15
tol
double tol
Definition mesh_smoothing.cpp:23
atom_test
int atom_test
Atom test.
Definition plastic.cpp:121

Variable Documentation

◆ help

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

Definition at line 26 of file project_roughness.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow