muMECH  1.0
tests_verified.cpp
Go to the documentation of this file.
1 #include "problem.h"
2 #include "mesh.h"
3 #include "homogenization.h"
4 #include "comparison.h"
5 #include "types.h"
6 #include "melement.h"
7 #include "vtk.h"
8 #include "tests.h"
9 #include "librw.h"
10 
11 
12 using namespace mumech;
13 
14 // ******************************************** NOVE ********************************************
16 int x3D_1I_Ellipsoid (void);
18 int x3D_1I_Ellipsoid_direct_API (void);
20 int x2D_1I_Ellipse_direct_API (void);
21 
23 int plstrain_2I_2d (void);
25 int plstrain_2I_3d (void);
26 
28 int grid_2d_3x3_3d (void);
29 int grid_2d_3x3_3d_F (void);
30 
32 int Four_3d_inclusions (void);
33 
34 // old
35 int Ellipsoid (void);
36 int Ellipse (void);
37 
38 int X2D_3incl (void);
39 int X3D_3incl (void);
40 
41 int x2D_1I_Ellipse_homog (void);
42 int x3D_1I_Ellipsoid_homog (void);
43 
44 //int integral_fields (void);
45 
46 
47 
49 void tests_verified (void)
50 {
51  int status = 0;
52 
53 
54  if (1) {
55  status += x3D_1I_Ellipsoid ();
56  status += x3D_1I_Ellipsoid_direct_API ();
57  status += x2D_1I_Ellipse_direct_API();
58 
59  status += plstrain_2I_2d ();
60  status += plstrain_2I_3d ();
61 
62  status += Four_3d_inclusions();
63 
64  status += grid_2d_3x3_3d ();
65  status += grid_2d_3x3_3d_F ();
66 
67  status += Ellipsoid ();
68  status += Ellipse ();
69 
70  status += X2D_3incl();
71  status += X3D_3incl();
72 
73  status += x2D_1I_Ellipse_homog();
74  status += x3D_1I_Ellipsoid_homog();
75 
76  //status += integral_fields ();
77  }
78 
79  if (status > 0) printf( "\n ERROR %d examples differs\n\n", status);
80  if (status == 0) printf( "\n OK all examples are ok\n\n");
81 }
82 
83 
85 int compare_fields (int n, const double *f1, const double *f2, double eps)
86 {
87  int nd = 0;
88 
89  while (n--)
90  if (fabs( (f1[n] - f2[n]) / f1[n]) > eps)
91  nd++;
92 
93  return nd;
94 }
95 
96 //
97 int diff_fields (int numDisp, int numStrain, int numStress, const double *d1, const double *d2, const double *e1, const double *e2, const double *s1, const double *s2, double eps)
98 {
99  int nd = 0;
100 
101  nd += compare_fields (numDisp, d1, d2, eps);
102  nd += compare_fields (numStrain, e1, e2, eps);
103  nd += compare_fields (numStress, s1, s2, eps);
104 
105  //if (nd) _warningg("error");
106  return nd;
107 }
108 
109 //
110 int cmp_vtk_file_delete_first (const char *file1, const char *file2)
111 {
112  bool status;
113 
114  if (file2)
115  status = cmp_vtk_file (file1, file2, 2);
116  else {
117  int l = strlen(file1);
118  char *file3 = new char[l+4];
119  strcpy(file3, file1);
120  sprintf(file3+l-4, ".OK");
121  strcpy(file3+l-1, file1+l-4);
122 
123  status = cmp_vtk_file (file1, file3, 2);
124  delete [] file3;
125  }
126 
127  // delete first file if files are same
128  if (status == false)
129  file_delete_if_exist_exit (file1);
130 
131  return (status ? 1 : 0);
132 }
133 
134 //
135 int cmp_txt_file_delete_first (const char *file1, const char *file2)
136 {
137  bool status;
138 
139  if (file2)
140  status = FP_cmp_files(file1, file2);
141  else {
142  int l = strlen(file1);
143  char *file3 = new char[l+4];
144  strcpy(file3, file1);
145  sprintf(file3+l-4, ".OK");
146  strcpy(file3+l-1, file1+l-4);
147 
148  status = FP_cmp_files(file1, file3);
149  delete [] file3;
150  }
151 
152  // delete first file if files are same
153  if (status == false)
154  file_delete_if_exist_exit (file1);
155 
156  return (status ? 1 : 0);
157 }
158 
159 
160 // **************************************************************************************************
161 // ******************************************** NOVE ********************************************
162 // **************************************************************************************************
163 
165 int x2D_1I_Ellipse_direct_API (void)
166 {
167  int status = 0;
168 
169  char vtkTopoFile[] = "verified/x2D_1I_Ellipse.topo.vtk";
170  //char vtkAmndFile[] = "verified/x2D_1I_Ellipse.amnd.vtk";
171  //char vtkAmndFileOK[] = "verified/x2D_1I_Ellipse.amnd.OK.vtk";
172  char vtkMeshFileOutPert[] = "verified/x2D_1I_Ellipse.rslt.pert";
173  char vtkMeshFileOutTotl[] = "verified/x2D_1I_Ellipse.rslt.totl";
174 
175  char vtkRsltFile[255];
176 
177  Problem *p = new Problem;
178 
179  // ************** DATA INPUT VIA FUNCTIONS ******************
180  // set problem
181  p->set_data_set();
182  p->set_dimension(2);
183  p->set_diffType (DT_NUMERICAL);
184  p->set_matrix_E_nu(1.0, 0.4);
185  //
186  p->set_numberOfRemoteStrains(4);
187  double rs1[] = { 1.0, 0.0, 0.0, 0.0 }; p->set_RemoteStrain(0, rs1);
188  double rs2[] = { 0.0, 0.0, 0.0, 1.0 }; p->set_RemoteStrain(1, rs2);
189  double rs3[] = { 0.0, 1.0, 1.0, 0.0 }; p->set_RemoteStrain(2, rs3);
190  double rs4[] = { 1.0, 3.0, 3.0, 2.0 }; p->set_RemoteStrain(3, rs4);
191 
192  p->set_number_of_inclusions(1);
193  // 1. inclusion
194  p->set_inclusion_centroids(0, 1.0, -2.0);
195  p->set_inclusion_E_nu (0, 2.4, 0.3);
196  p->set_inclusion_semiaxesDimensions(0, 2.0, 1.0);
197  p->set_inclusion_EulerAnglesDEG(0, 70.0);
198 
199  p->input_data_initialize_and_check_consistency();
200  // ************** DATA INPUT VIA FUNCTIONS ******************
201 
202  p->convert_to_equivalent_problem ();
203  //p->print_equivalent_problem (vtkAmndFile);
204  p->print_visualization(vtkTopoFile, 20, 2);
205 
206  // regular mesh, cube 8x8x8 with centroid 1.0 -2.0 4.0
207  double p1[3] = { 5.0, 2.0, 0.0};
208  double p2[3] = {-3.0, -6.0, 0.0};
209  long n[3] = { 8, 8, 0};
210 
211  //status += cmp_vtk_file_delete_first (vtkAmndFile, vtkAmndFileOK);
212 
213  // PERTURBATION
214  p->printFieldsOnMeshGrid(vtkMeshFileOutPert, p1, p2, n, 'p', 0, 0, PFCM_FULL);
215 
216  // compare and delete files
217  for (int i=0; i<p->give_nLC(); i++) {
218  sprintf(vtkRsltFile, "%s.%02d.vtk", vtkMeshFileOutPert, i);
219  status += cmp_vtk_file_delete_first (vtkRsltFile);
220  }
221 
222  // TOTAL
223  p->printFieldsOnMeshGrid(vtkMeshFileOutTotl, p1, p2, n, 't', 0, 0, PFCM_FULL);
224 
225  // compare and delete files
226  for (int i=0; i<p->give_nLC(); i++) {
227  sprintf(vtkRsltFile, "%s.%02d.vtk", vtkMeshFileOutTotl, i);
228  status += cmp_vtk_file_delete_first (vtkRsltFile);
229  }
230 
231  //
232  delete p;
233 
234  return status;
235 }
236 
237 
239 int x3D_1I_Ellipsoid_direct_API (void)
240 {
241  int status = 0;
242  char vtkAmndFile[] = "verified/x3D_1I_Ellipsoid.amnd.vtk";
243  char vtkMeshFileOutPert[] = "verified/x3D_1I_Ellipsoid.rslt.pert";
244  char vtkMeshFileOutTotl[] = "verified/x3D_1I_Ellipsoid.rslt.totl";
245 
246  char vtkRsltFile[255];
247 
248  Problem *p = new Problem;
249 
250  // ************** DATA INPUT ******************
251  // START - DATA INPUT VIA FUNCTIONS
252  // do not change the order of command
253 
254  // set problem
255  p->set_data_set();
256  p->set_dimension(3);
257  p->set_diffType (DT_ANALITICAL);
258  p->set_matrix_E_nu(1.0, 0.4);
259  //
260  p->set_numberOfRemoteStrains(7);
261  double rs1[] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; p->set_RemoteStrain(0, rs1);
262  double rs2[] = { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 }; p->set_RemoteStrain(1, rs2);
263  double rs3[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0 }; p->set_RemoteStrain(2, rs3);
264  double rs4[] = { 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; p->set_RemoteStrain(3, rs4);
265  double rs5[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0 }; p->set_RemoteStrain(4, rs5);
266  double rs6[] = { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0 }; p->set_RemoteStrain(5, rs6);
267  double rs7[] = { 1.0, 4.0, 6.0, 4.0, 2.0, 5.0, 6.0, 5.0, 3.0 }; p->set_RemoteStrain(6, rs7);
268 
269  // loop over all inclusions
270  p->set_number_of_inclusions(1);
271  // 1. inclusion
272  p->set_inclusion_centroids(0,1.0, -2.0, 4.0);
273  p->set_inclusion_E_nu(0, 2.4, 0.3);
274  p->set_inclusion_semiaxesDimensions(0, 2.0, 1.0, 0.5);
275  p->set_inclusion_EulerAnglesDEG(0, 20.0, 40.0, 50.0);
276 
277  p->input_data_initialize_and_check_consistency();
278  // ************** DATA INPUT ******************
279 
280  p->convert_to_equivalent_problem ();
281 
282  p->print_equivalent_problem (vtkAmndFile);
283 
284 
285  // regular mesh, cube 8x8x8 with centroid 1.0 -2.0 4.0
286  double p1[3] = { 5.0, 2.0, 8.0};
287  double p2[3] = {-3.0, -6.0, 0.0};
288  long n[3] = { 8, 8, 8};
289 
290  status += cmp_vtk_file_delete_first (vtkAmndFile);
291 
292  // PERTURBATION
293  p->printFieldsOnMeshGrid(vtkMeshFileOutPert, p1, p2, n, 'p', 0, 0, PFCM_FULL);
294 
295  // compare and delete files
296  for (int i=0; i<p->give_nLC(); i++) {
297  sprintf(vtkRsltFile, "%s.%02d.vtk", vtkMeshFileOutPert, i);
298  status += cmp_vtk_file_delete_first (vtkRsltFile);
299  }
300 
301  // TOTAL
302  p->printFieldsOnMeshGrid(vtkMeshFileOutTotl, p1, p2, n, 't', 0, 0, PFCM_FULL);
303 
304  // compare and delete files
305  for (int i=0; i<p->give_nLC(); i++) {
306  sprintf(vtkRsltFile, "%s.%02d.vtk", vtkMeshFileOutTotl, i);
307  status += cmp_vtk_file_delete_first (vtkRsltFile);
308  }
309 
310  //
311  delete p;
312  return status;
313 }
314 
315 
316 
318 int x3D_1I_Ellipsoid (void)
319 {
320  int status = 0;
321  char vtkGeomFile[] = "verified/x3D_1I_Ellipsoid.geom.vtk";
322  char vtkAmndFile[] = "verified/x3D_1I_Ellipsoid.amnd.vtk";
323  char vtkMeshFileOut[] = "verified/x3D_1I_Ellipsoid.rslt.pert";
324 
325  char vtkRsltFile[255];
326 
327  Problem *p = new Problem;
328 
329  // regular mesh, cube 8x8x8 with centroid 1.0 -2.0 4.0
330  double p1[3] = { 5.0, 2.0, 8.0};
331  double p2[3] = {-3.0, -6.0, 0.0};
332  long n[3] = { 8, 8, 8};
333 
334  // compute
335  p->read_inclusions_plus_initialize_and_print (vtkGeomFile, vtkAmndFile, DT_ANALITICAL);
336  p->printFieldsOnMeshGrid(vtkMeshFileOut, p1, p2, n, 'p', 0, 0, PFCM_FULL);
337 
338  // compare files
339  status += cmp_vtk_file_delete_first (vtkAmndFile);
340  for (int i=0; i<p->give_nLC(); i++) {
341  sprintf(vtkRsltFile, "%s.%02d.vtk", vtkMeshFileOut, i);
342  status += cmp_vtk_file_delete_first (vtkRsltFile);
343  }
344 
345  //
346  delete p;
347  return status;
348 }
349 
351 int grid_2d_3x3_3d (void)
352 {
353  int status = 0;
354  char vtkAmndFile[] = "verified/grid_2d_3x3_3d.amnd.vtk";
355  char vtkTopoFile[] = "verified/grid_2d_3x3_3d.topo.vtk";
356  char vtkMeshFile[] = "verified/grid_2d_3x3_3d.rslt";
357  char vtkMeshFileOut[] = "verified/grid_2d_3x3_3d.rslt.00.vtk";
358 
359  Problem *p = new Problem;
360 
361  // ************** DATA INPUT VIA FUNCTIONS ******************
362  // problem
363  p->set_DataDimDiff (3, DT_ANALITICAL);
364  p->set_SBA_optimized(true);
365  //
366  double rs [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
367  p->set_matrix_E_nu(1.0, 0.3);
368  p->set_numberOfRemoteStrains(1);
369  p->set_RemoteStrain(0, rs);
370  // inclusions
371  p->set_number_of_inclusions(9);
372  // inclID x y z E nu a1 a2 a3 ae1 ae2 ae3
373  p->set_inclusion_all (0, 0.0, 0.0, 0.0, 10, 0.4, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0);
374  p->set_inclusion_all (1, 4.0, 0.0, 0.0, 10, 0.4, 1.0, 0.8, 0.5, 0.0, 0.0, 0.0);
375  p->set_inclusion_all (2, 8.0, 0.0, 0.0, 10, 0.4, 1.0, 0.7, 0.4, 30.0, 40.0, 50.0);
376  p->set_inclusion_all (3, 0.0, 4.0, 0.0, 10, 0.4, 1.0, 1.0, 0.5, 0.0, 70.0, 0.0);
377  p->set_inclusion_all (4, 4.0, 4.0, 0.0, 10, 0.4, 1.0, 1.0, 0.6, 70.0, 40.0, 0.0);
378  p->set_inclusion_all (5, 8.0, 4.0, 0.0, 10, 0.4, 1.0, 0.5, 0.5, 0.0, 0.0, 90.0);
379  p->set_inclusion_all (6, 0.0, 8.0, 0.0, 10, 0.4, 1.0, 0.5, 0.3, 30.0, 30.0, 30.0);
380  p->set_inclusion_all (7, 4.0, 8.0, 0.0, 10, 0.4, 1.0, 1.0, 0.3, 0.0, 0.0, 0.0);
381  p->set_inclusion_all (8, 8.0, 8.0, 0.0, 10, 0.4, 1.0, 0.8, 0.5, 90.0, 0.0, 0.0);
382  //
383  p->input_data_initialize_and_check_consistency();
384  // ************** DATA INPUT VIA FUNCTIONS ******************
385  p->set_SBAM(SBAT_ORIGINAL);
386  p->convert_to_equivalent_problem();
387  p->print_equivalent_problem(vtkAmndFile);
388  p->print_visualization(vtkTopoFile, 7, 3);
389 
390  // regular mesh
391  double p1[3] = {-4.0, -4.0, 0.0 };
392  double p2[3] = {12.0, 12.0, 0.0 };
393  long nn[3] = {100, 100, 0 };
394 
395  p->set_intFieldsShape(4);
396  p->printFieldsOnMeshGrid(vtkMeshFile, p1, p2, nn, 'p', 0, 0, PFCM_FULL);
397 
398  status += cmp_vtk_file_delete_first (vtkAmndFile);
399  status += cmp_vtk_file_delete_first (vtkMeshFileOut);
400 
401  delete p;
402  return status;
403 }
405 int grid_2d_3x3_3d_F (void)
406 {
407  int status = 0;
408  char vtkAmndFile[] = "verified/grid_2d_3x3_3dF.amnd.vtk";
409  char vtkMeshFile[] = "verified/grid_2d_3x3_3dF.rslt";
410  char vtkMeshFileOut[] = "verified/grid_2d_3x3_3dF.rslt.00.vtk";
411 
412  Problem *p = new Problem;
413 
414  // ************** DATA INPUT VIA FUNCTIONS ******************
415  // problem
416  p->set_DataDimDiff (3, DT_ANALITICAL);
417  p->set_SBA_optimized(false);
418  //
419  double rs [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
420  p->set_matrix_E_nu(1.0, 0.3);
421  p->set_numberOfRemoteStrains(1);
422  p->set_RemoteStrain(0,rs);
423  // inclusions
424  p->set_number_of_inclusions(9);
425  // inclID x y z E nu a1 a2 a3 ae1 ae2 ae3
426  p->set_inclusion_all (0, 0.0, 0.0, 0.0, 10, 0.4, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0);
427  p->set_inclusion_all (1, 4.0, 0.0, 0.0, 10, 0.4, 1.0, 0.8, 0.5, 0.0, 0.0, 0.0);
428  p->set_inclusion_all (2, 8.0, 0.0, 0.0, 10, 0.4, 1.0, 0.7, 0.4, 30.0, 40.0, 50.0);
429  p->set_inclusion_all (3, 0.0, 4.0, 0.0, 10, 0.4, 1.0, 1.0, 0.5, 0.0, 70.0, 0.0);
430  p->set_inclusion_all (4, 4.0, 4.0, 0.0, 10, 0.4, 1.0, 1.0, 0.6, 70.0, 40.0, 0.0);
431  p->set_inclusion_all (5, 8.0, 4.0, 0.0, 10, 0.4, 1.0, 0.5, 0.5, 0.0, 0.0, 90.0);
432  p->set_inclusion_all (6, 0.0, 8.0, 0.0, 10, 0.4, 1.0, 0.5, 0.3, 30.0, 30.0, 30.0);
433  p->set_inclusion_all (7, 4.0, 8.0, 0.0, 10, 0.4, 1.0, 1.0, 0.3, 0.0, 0.0, 0.0);
434  p->set_inclusion_all (8, 8.0, 8.0, 0.0, 10, 0.4, 1.0, 0.8, 0.5, 90.0, 0.0, 0.0);
435  //
436  p->input_data_initialize_and_check_consistency();
437  // ************** DATA INPUT VIA FUNCTIONS ******************
438 
439  p->set_SBAM(SBAT_ORIGINAL);
440  p->convert_to_equivalent_problem();
441  p->print_equivalent_problem(vtkAmndFile);
442 
443  // regular mesh
444  double p1[3] = {-4.0, -4.0, 0.0 };
445  double p2[3] = {12.0, 12.0, 0.0 };
446  long nn[3] = {100, 100, 0 };
447 
448  p->set_intFieldsShape(4);
449  p->printFieldsOnMeshGrid(vtkMeshFile, p1, p2, nn, 'p', 0, 0, PFCM_FULL);
450 
451  status += cmp_vtk_file_delete_first (vtkAmndFile);
452  status += cmp_vtk_file_delete_first (vtkMeshFileOut);
453 
454  delete p;
455  return status;
456 }
457 
459 int plstrain_2I_2d (void)
460 {
461  Problem *p = new Problem;
462 
463  // ************** DATA INPUT VIA FUNCTIONS ******************
464  // set problem
465  p->set_data_set();
466  p->set_dimension(2);
467  p->set_diffType(DT_NUMERICAL);
468  p->set_matrix_E_nu(1.0, 0.2);
469  //
470  p->set_numberOfRemoteStrains(1);
471  double rs [] = { 1.0, 0.0, 0.0, 0.0 };
472  p->set_RemoteStrain(0, rs);
473 
474  p->set_number_of_inclusions(2);
475  // 1. inclusion
476  p->set_inclusion_centroids(0,-2.0, 0.0);
477  p->set_inclusion_E_nu(0, 10, 0.4);
478  p->set_inclusion_semiaxesDimensions(0, 1.0, 0.5);
479  p->set_inclusion_EulerAnglesDEG(0, 0.0);
480  // 2. inclusion
481  p->set_inclusion_centroids(1,2.0, 0.0);
482  p->set_inclusion_E_nu(1, 10, 0.4);
483  p->set_inclusion_semiaxesDimensions(1, 1.0, 0.5);
484  p->set_inclusion_EulerAnglesDEG(1, 0.0);
485 
486  p->input_data_initialize_and_check_consistency();
487  // ************** DATA INPUT VIA FUNCTIONS ******************
488 
489  //p->set_SBA_requiredNumberOfIterations ( 0); p->set_addExt(false);
490 
491  // Arbitrary coordinates of one point
492  double coords1[3] = { 2.0, 0.0, 0.0}; // inside point
493  double coords2[3] = { 0.0, 0.0, 0.0}; // outside point
494 
495  // Allocate arrays
496  const int numDisp = 2;
497  const int numStrain = 4;
498  const int numStress = 4;
499  // Reference calculated fields
500  double d1ok[numDisp], d1[numDisp], d2ok[numDisp], d2[numDisp];
501  double e1ok[numStrain], e1[numStrain], e2ok[numStrain], e2[numStrain];
502  double s1ok[numStrain], s1[numStrain], s2ok[numStrain], s2[numStrain];
503 
504  char vtkAmndFile [] = "verified/plstrain_2I_2d.vtk";
505 
506  p->set_SBA_optimized(false);
507  p->convert_to_equivalent_problem();
508  p->print_equivalent_problem (vtkAmndFile);
509  p->giveFieldsOfPointOneRS(d1ok, e1ok, s1ok, coords1, 'p', 0, PFCM_OPTIMIZED);
510  p->giveFieldsOfPointOneRS(d2ok, e2ok, s2ok, coords2, 'p', 0, PFCM_OPTIMIZED);
511  delete p;
512 
513  p = new Problem;
514  p->read_input_file (vtkAmndFile);
515  p->input_data_initialize_and_check_consistency();
516  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_OPTIMIZED);
517  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_OPTIMIZED);
518  delete p;
519 
520  int status = 0;
521  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) {
522  _warningg("2 ellipses results differs"); status++;
523  }
524  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) {
525  _warningg("2 ellipses results differs"); status++;
526  }
527  status += cmp_vtk_file_delete_first (vtkAmndFile);
528  return status;
529 }
530 
532 int plstrain_2I_3d (void)
533 {
534  Problem *p = new Problem;
535 
536  // ************** DATA INPUT VIA FUNCTIONS ******************
537  // set problem
538  p->set_data_set();
539  p->set_dimension(3);
540  p->set_diffType(DT_ANALITICAL);
541  p->set_matrix_E_nu(1.0, 0.2);
542  //
543  p->set_numberOfRemoteStrains(1);
544  double rs [] = { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
545  p->set_RemoteStrain(0, rs);
546 
547  p->set_number_of_inclusions(2);
548  // 1. inclusion
549  p->set_inclusion_centroids(0,-2.0, 0.0, 0.0);
550  p->set_inclusion_E_nu(0, 10, 0.4);
551  p->set_inclusion_semiaxesDimensions(0, 100*1.0, 1.0, 0.5);
552  p->set_inclusion_EulerAnglesDEG(0, 0.0, 90.0, 90.0);
553  // 2. inclusion
554  p->set_inclusion_centroids(1,2.0, 0.0, 0.0);
555  p->set_inclusion_E_nu(1, 10, 0.4);
556  p->set_inclusion_semiaxesDimensions(1, 100*1.0, 1.0, 0.5);
557  p->set_inclusion_EulerAnglesDEG(1, 0.0, 90.0, 90.0);
558 
559  p->input_data_initialize_and_check_consistency();
560  // ************** DATA INPUT VIA FUNCTIONS ******************
561 
562  //p->set_SBA_requiredNumberOfIterations ( 0); p->set_addExt(false);
563 
564  // Arbitrary coordinates of one point
565  double coords1[3] = { 2.0, 0.0, 0.0}; // inside point
566  double coords2[3] = { 0.0, 0.0, 0.0}; // outside point
567 
568  // Allocate arrays
569  const int numDisp = 3;
570  const int numStrain = 9;
571  const int numStress = 9;
572  // Reference calculated fields
573  double d1ok[numDisp], d1[numDisp], d2ok[numDisp], d2[numDisp];
574  double e1ok[numStrain], e1[numStrain], e2ok[numStrain], e2[numStrain];
575  double s1ok[numStrain], s1[numStrain], s2ok[numStrain], s2[numStrain];
576 
577  char vtkAmndFile [] = "verified/plstrain_2I_3d.vtk";
578 
579  p->set_SBA_optimized(false);
580  p->convert_to_equivalent_problem();
581  p->print_equivalent_problem (vtkAmndFile);
582  p->giveFieldsOfPointOneRS(d1ok, e1ok, s1ok, coords1, 'p', 0, PFCM_OPTIMIZED);
583  p->giveFieldsOfPointOneRS(d2ok, e2ok, s2ok, coords2, 'p', 0, PFCM_OPTIMIZED);
584  delete p;
585 
586  p = new Problem;
587  p->read_input_file (vtkAmndFile);
588  p->input_data_initialize_and_check_consistency();
589  p->set_diffType(DT_ANALITICAL);
590  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_OPTIMIZED);
591  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_OPTIMIZED);
592  delete p;
593 
594  int status = 0;
595  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("2 ellipses results differs"); status++; }
596  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("2 ellipses results differs"); status++; }
597  status += cmp_vtk_file_delete_first (vtkAmndFile);
598  return status;
599 }
600 
602 
604 int x2D_1I_Ellipse_homog(void)
605 {
606  int status = 0;
607  char rsltsFile[] = "verified/homogenization_1I_2d.txt";
608 
609  Problem *p = new Problem;
610 
611  // ************** DATA INPUT VIA FUNCTIONS ******************
612  // set problem
613  p->set_DataDimDiff(2, DT_NUMERICAL);
614  p->set_matrix_E_nu (10.0, 0.3);
615  p->set_UnitRemoteStrains();
616 
617  p->set_number_of_inclusions(1);
618  // 1. inclusion E nu
619  p->set_inclusion_all(0, 0.0, 0.0, 1.0, 0.2, 2.0, 1.0, 0.0);
620  // ************** DATA INPUT VIA FUNCTIONS ******************
621 
622  p->input_data_initialize_and_check_consistency();
623  p->convert_to_equivalent_problem ();
624 
625  Homogenization *h;
626 
627  double sm[3*3];
628  double x = 2.3;
629  double y = 1.3;
630 
631  FILE *file = fopen(rsltsFile, "w");
632 
633  // *** STANDARDNI HOMOGENIZACE ***
634  p->matrix_giveFullStiffMatrix(sm); fprintf (file, "%s:\t\t", "Matrix" ); for (int i=0; i<9; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
635  h = p->give_new_homogenization (HT_VoightBound); h->set_boundingBox (-x, -y, x, y); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t", h->giveClassName()); for (int i=0; i<9; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
636  h = p->give_new_homogenization (HT_RegGrid); h->set_boundingBox (-x, -y, x, y); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t", h->giveClassName()); for (int i=0; i<9; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
637 
638  // not implemented
639  //h = p->give_new_homogenization (HT_SelfConsistent); h->set_boundingBox (-x, -y, x, y); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t", h->giveClassName()); for (int i=0; i<9; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
640  // 3d only
641  //h = p->give_new_homogenization (HT_DifferentialScheme); h->set_boundingBox (-x, -y, x, y); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t", h->giveClassName()); for (int i=0; i<9; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
642 
643  h = p->give_new_homogenization (HT_MoriTanaka); h->set_boundingBox (-x, -y, x, y); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t", h->giveClassName()); for (int i=0; i<9; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
644  h = p->give_new_homogenization (HT_Dilute); h->set_boundingBox (-x, -y, x, y); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t\t", h->giveClassName()); for (int i=0; i<9; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
645  h = p->give_new_homogenization (HT_ReussBound); h->set_boundingBox (-x, -y, x, y); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t", h->giveClassName()); for (int i=0; i<9; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
646  p->give_inclusion(0)->give_StiffnessMatrixFull(sm); fprintf (file, "%s:\t", "Inclusion" ); for (int i=0; i<9; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
647 
648  double fv = h->giveTotalVolumeFractionOfInclusions();
649  fprintf (file, "Fraction volume: %6.2lf\n", 100*fv);
650 
651  fclose(file);
652 
653  delete p;
654 
655  // compare results
656  status += cmp_txt_file_delete_first (rsltsFile);
657 
658  return status;
659 }
660 
662 int x3D_1I_Ellipsoid_homog(void)
663 {
664  int status = 0;
665  char rsltsFile[] = "verified/homogenization_1I_3d.txt";
666 
667  Problem *p = new Problem;
668 
669  // ************** DATA INPUT VIA FUNCTIONS ******************
670  // set problem
671  p->set_DataDimDiff(3, DT_ANALITICAL);
672  p->set_matrix_E_nu (10.0, 0.2);
673  p->set_UnitRemoteStrains();
674 
675  p->set_number_of_inclusions(1);
676  // 1. inclusion E nu
677  p->set_inclusion_all(0, 0.0, 0.0, 0.0, 1.0, 0.2, 4.0, 2.0, 1.0, 0.0, 90.0, 90.0);
678  // ************** DATA INPUT VIA FUNCTIONS ******************
679 
680  p->input_data_initialize_and_check_consistency();
681  p->convert_to_equivalent_problem ();
682 
683  Homogenization *h;
684 
685  double sm[6*6];
686  double x = 2.5;
687  double y = 1.5;
688  double z = 5.0;
689 
690  FILE *file = fopen(rsltsFile, "w");
691 
692  // *** STANDARDNI HOMOGENIZACE ***
693  p->matrix_giveFullStiffMatrix(sm); fprintf (file, "%s:\t\t\t", "Matrix" ); for (int i=0; i<36; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
694  h = p->give_new_homogenization (HT_VoightBound); h->set_boundingBox (-x, -y, -z, x, y, z); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t\t", h->giveClassName()); for (int i=0; i<36; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
695  h = p->give_new_homogenization (HT_RegGrid); h->set_boundingBox (-x, -y, -z, x, y, z); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t\t", h->giveClassName()); for (int i=0; i<36; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
696  h = p->give_new_homogenization (HT_DifferentialScheme); h->set_boundingBox (-x, -y, -z, x, y, z); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t", h->giveClassName()); for (int i=0; i<36; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
697  h = p->give_new_homogenization (HT_MoriTanaka); h->set_boundingBox (-x, -y, -z, x, y, z); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t\t", h->giveClassName()); for (int i=0; i<36; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
698  h = p->give_new_homogenization (HT_Dilute); h->set_boundingBox (-x, -y, -z, x, y, z); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t\t\t", h->giveClassName()); for (int i=0; i<36; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
699  h = p->give_new_homogenization (HT_ReussBound); h->set_boundingBox (-x, -y, -z, x, y, z); h->giveHomogenizedStiffnessMatrix (sm); fprintf (file, "%s:\t\t", h->giveClassName()); for (int i=0; i<36; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
700  p->give_inclusion(0)->give_StiffnessMatrixFull(sm); fprintf (file, "%s:\t\t", "Inclusion" ); for (int i=0; i<36; i++) fprintf (file, " %15.8lf", sm[i]); fprintf (file, "\n");
701 
702  double fv = h->giveTotalVolumeFractionOfInclusions();
703  fprintf (file, "Fraction volume: %6.2lf\n", 100*fv);
704 
705  fclose(file);
706 
707  delete p;
708 
709  // compare results
710  status += cmp_txt_file_delete_first (rsltsFile);
711 
712  return status;
713 }
714 
715 
716 
717 
718 
719 
721 // ************************** STARE *************************************************************
723 
724 
725 
726 // Four 3d inclusions.
727 int Four_3d_inclusions (void)
728 {
729  int status = 0;
730  char vtkGeomFile[] = "verified/four_3d_icls.geom.vtk";
731  char vtkAmndFile[] = "verified/four_3d_icls.amnd.vtk";
732  char vtkTopoFile[] = "verified/four_3d_icls.topo.vtk";
733  char vtkRsltFileA[] = "verified/four_3d_icls.rsltA";
734  char vtkRsltFileD[] = "verified/four_3d_icls.rsltD";
735  char vtkRsltFilevtkA[] = "verified/four_3d_icls.rsltA.00.vtk";
736  char vtkRsltFilevtkD[] = "verified/four_3d_icls.rsltD.00.vtk";
737 
738  // regular mesh
739  double p1[3] = { -5.0, -8.0, -5.0 };
740  double p2[3] = { 10.0, 8.0, 9.0 };
741  long nn[3] = { 30, 30, 30};
742 
743  Problem *p;
744 
745  // CREATE AMENDED FILE + COMPUTE FIELDS ON MESH
746  p = new Problem;
747  p->set_diffType (DT_ANALITICAL);
748  p->read_input_file (vtkGeomFile);
749  p->input_data_initialize_and_check_consistency();
750  p->set_SBAM(SBAT_ORIGINAL);
751  p->convert_to_equivalent_problem ();
752  p->print_equivalent_problem (vtkAmndFile);
753  delete p;
754 
755  p = new Problem;
756  p->set_diffType (DT_ANALITICAL);
757  p->read_input_file (vtkAmndFile);
758  p->input_data_initialize_and_check_consistency();
759  p->set_intFieldsShape(4);
760  p->printFieldsOnMeshGrid(vtkRsltFileA, p1, p2, nn, 'p', 0, 0, PFCM_FULL);
761  delete p;
762 
763  // compare files
764  status += cmp_vtk_file_delete_first (vtkAmndFile);
765  status += cmp_vtk_file_delete_first (vtkRsltFilevtkA);
766 
767  // DIRECT COMPUTATION FIELDS ON MESH (+ PRINT AMENDED FILE)
768  p = new Problem;
769  p->set_SBAM(SBAT_ORIGINAL);
770  p->read_inclusions_plus_initialize_and_print (vtkGeomFile, NULL, DT_ANALITICAL);
771  p->print_visualization(vtkTopoFile, 10, 3);
772  p->print_equivalent_problem (vtkAmndFile);
773  p->set_intFieldsShape(4);
774  p->printFieldsOnMeshGrid(vtkRsltFileD, p1, p2, nn, 'p', 0, 0, PFCM_FULL);
775  delete p;
776 
777  // compare files
778  status += cmp_vtk_file_delete_first (vtkAmndFile);
779  status += cmp_vtk_file_delete_first (vtkRsltFilevtkD);
780  return status;
781 }
782 
783 //
784 int Ellipse (void)
785 {
786  int status = 0;
787  // Inclusion's geometry and topology files
788  char vtkGeomFile[] = "verified/Ellipse.geom.vtk";
789  char vtkGeomFileField[] = "verified/Ellipse.geom.field.vtk";
790  //char vtuGeomFile[] = "verified/Ellipsoid.geom.vtu";
791  //char vtuGeomFileField[] = "verified/Ellipsoid.geom.field.vtu";
792  char vtkAmndFile[] = "verified/Ellipse.amnd.vtk";
793  //char vtkAmndFileOK[] = "verified/Ellipse.amnd.OK.vtk";
794 
795  // Infinite medium properties
796  double E = 1.0; // Young's modulus of matrix
797  double nu = 0.4; // Poisson number of matrix
798  bool sbaopt = true; // type of self-balancing algorithm _OPTIMIZED_
799 
800  // Arbitrary coordinates of one point
801  double coords1[3] = { 1.4, 2.4, 0.0}; // inside point
802  double coords2[3] = {-1.1, 4.1, 0.0}; // outside point
803 
804  Problem *p;
805 
806  // Allocate arrays
807  const int numDisp = 2;
808  const int numStrain = 4;
809  const int numStress = 4;
810  // Reference calculated fields
811  double d1ok[numDisp], d1[numDisp], d2ok[numDisp], d2[numDisp];
812  double e1ok[numStrain], e1[numStrain], e2ok[numStrain], e2[numStrain];
813  double s1ok[numStrain], s1[numStrain], s2ok[numStrain], s2[numStrain];
814 
815  // *******************************
816  // *** *** *** VTK *** *** ***
817  // *******************************
818  // nactu core file, udelam z nej rec; nactu rec a spocitam
819  // sbalg a matrix se zadavaji rucne
820  if (1) {
821  p = new Problem;
822  p->read_input_file (vtkGeomFile);
823  p->set_SBA_optimized(sbaopt);
824  p->set_matrix_E_nu(E, nu);
825  p->input_data_initialize_and_check_consistency();
826  p->convert_to_equivalent_problem ();
827  p->print_equivalent_problem (vtkAmndFile);
828  p->giveFieldsOfPointOneRS(d1ok, e1ok, s1ok, coords1, 'p', 0, PFCM_OPTIMIZED);
829  p->giveFieldsOfPointOneRS(d2ok, e2ok, s2ok, coords2, 'p', 0, PFCM_OPTIMIZED);
830  delete p;
831 
832  p = new Problem;
833  p->read_input_file (vtkAmndFile);
834  p->input_data_initialize_and_check_consistency();
835  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_OPTIMIZED);
836  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_OPTIMIZED);
837  delete p;
838 
839  // termitovo: NEFUNGUJE REC FILE, PROTO JE NASLEDUJICI TEST ZAKOMENTOVANEJ !!!!!!!!!!!!!
840 
841  //if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
842  //if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
843 
844  // compare and delete files
845  //status += cmp_vtk_file_delete_first (vtkAmndFile, vtkAmndFileOK);
846  file_delete_if_exist(vtkAmndFile);
847  }
848  // nactu core file a spocitam
849  if (1) {
850  p = new Problem;
851  p->read_input_file (vtkGeomFileField);
852  p->input_data_initialize_and_check_consistency();
853  p->convert_to_equivalent_problem ();
854  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_OPTIMIZED);
855  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_OPTIMIZED);
856  delete p;
857 
858  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
859  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
860  }
861  return status;
862 }
863 
864 
865 //
866 int Ellipsoid (void)
867 {
868  int status = 0;
869  // Inclusion's geometry and topology files
870  char vtkGeomFile[] = "verified/Ellipsoid.geom.vtk";
871  char vtkGeomFileField[] = "verified/Ellipsoid.geom.field.vtk";
872  char vtuGeomFile[] = "verified/Ellipsoid.geom.vtu";
873  char vtuGeomFileField[] = "verified/Ellipsoid.geom.field.vtu";
874  char vtkAmndFile[] = "verified/Ellipsoid.amnd.vtk";
875 
876  // Infinite medium properties
877  double E = 1.0; // Young's modulus of matrix
878  double nu = 0.4; // Poisson number of matrix
879  bool sbaopt = true; // type of self-balancing algorithm _OPTIMIZED_
880 
881  // Arbitrary coordinates of one point
882  double coords1[3] = { 1.4, 2.4, 5.4}; // inside point
883  double coords2[3] = {-1.1, 4.1, 2.9}; // outside point
884 
885  Problem *p;
886 
887  // Allocate arrays
888  const int numDisp = 3;
889  const int numStrain = 9;
890  const int numStress = 9;
891  // Reference calculated fields
892  double d1ok[numDisp], d1[numDisp], d2ok[numDisp], d2[numDisp];
893  double e1ok[numStrain], e1[numStrain], e2ok[numStrain], e2[numStrain];
894  double s1ok[numStrain], s1[numStrain], s2ok[numStrain], s2[numStrain];
895 
896  // *******************************
897  // *** *** *** VTK *** *** ***
898  // *******************************
899  // nactu core file, udelam z nej rec; nactu rec a spocitam
900  // sbalg a matrix se zadavaji rucne
901  if (1) {
902  p = new Problem;
903  p->read_input_file (vtkGeomFile);
904  p->set_SBA_optimized(sbaopt);
905  p->set_matrix_E_nu(E, nu);
906  p->input_data_initialize_and_check_consistency();
907  p->convert_to_equivalent_problem ();
908  p->print_equivalent_problem (vtkAmndFile);
909  delete p;
910 
911  p = new Problem;
912  p->read_input_file (vtkAmndFile);
913  p->input_data_initialize_and_check_consistency();
914  p->set_diffType (DT_ANALITICAL);
915  p->giveFieldsOfPointOneRS(d1ok, e1ok, s1ok, coords1, 'p', 0, PFCM_FULL);
916  p->giveFieldsOfPointOneRS(d2ok, e2ok, s2ok, coords2, 'p', 0, PFCM_FULL);
917  delete p;
918 
919  // compare and delete files
920  status += cmp_vtk_file_delete_first (vtkAmndFile);
921  }
922  // nactu core file, udelam z nej rec; nactu rec a spocitam
923  // sbalg a matrix se zadavaji v souboru
924  if (1) {
925  p = new Problem;
926  p->read_input_file (vtkGeomFileField);
927  p->input_data_initialize_and_check_consistency();
928  p->convert_to_equivalent_problem ();
929  p->print_equivalent_problem (vtkAmndFile);
930  delete p;
931 
932  p = new Problem;
933  p->read_input_file (vtkAmndFile);
934  p->input_data_initialize_and_check_consistency();
935  p->set_diffType (DT_ANALITICAL);
936  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_FULL);
937  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_FULL);
938  delete p;
939 
940  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
941  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
942  // compare and delete files
943  status += cmp_vtk_file_delete_first (vtkAmndFile);
944  }
945  // nactu core file a spocitam
946  if (1) {
947  p = new Problem;
948  p->read_input_file (vtkGeomFileField);
949  p->input_data_initialize_and_check_consistency();
950  p->convert_to_equivalent_problem ();
951  p->set_diffType (DT_ANALITICAL);
952  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_FULL);
953  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_FULL);
954  delete p;
955 
956  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
957  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
958  }
959  // nactu core file, vytisknu rec a spocitam
960  if (1) {
961  p = new Problem;
962  p->read_input_file (vtkGeomFileField);
963  p->input_data_initialize_and_check_consistency();
964  p->convert_to_equivalent_problem ();
965  p->print_equivalent_problem (vtkAmndFile);
966  p->set_diffType (DT_ANALITICAL);
967  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_FULL);
968  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_FULL);
969  delete p;
970 
971  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
972  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
973  // compare and delete files
974  status += cmp_vtk_file_delete_first (vtkAmndFile);
975  }
976  // nactu core file, spocitam a vytisknu rec
977  if (1) {
978  p = new Problem;
979  p->read_input_file (vtkGeomFileField);
980  p->input_data_initialize_and_check_consistency();
981  p->convert_to_equivalent_problem ();
982  p->set_diffType (DT_ANALITICAL);
983  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_FULL);
984  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_FULL);
985  p->print_equivalent_problem (vtkAmndFile);
986  delete p;
987 
988  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
989  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
990  // compare and delete files
991  status += cmp_vtk_file_delete_first (vtkAmndFile);
992  }
993 
994 
995  // *******************************
996  // *** *** *** VTU *** *** ***
997  // *******************************
998  // nactu core file, udelam z nej rec; nactu rec a spocitam
999  // sbalg a matrix se zadavaji rucne
1000  if (1) {
1001  p = new Problem;
1002  p->read_input_file (vtuGeomFile);
1003  p->set_SBA_optimized(sbaopt);
1004  p->set_matrix_E_nu(E, nu);
1005  p->input_data_initialize_and_check_consistency();
1006  p->convert_to_equivalent_problem ();
1007  p->print_equivalent_problem (vtkAmndFile);
1008  delete p;
1009 
1010  //
1011  p = new Problem;
1012  p->read_input_file (vtkAmndFile);
1013  p->input_data_initialize_and_check_consistency();
1014  p->set_diffType (DT_ANALITICAL);
1015  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_FULL);
1016  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_FULL);
1017  delete p;
1018 
1019  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1020  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1021  // compare and delete files
1022  status += cmp_vtk_file_delete_first (vtkAmndFile);
1023  }
1024  // nactu core file, udelam z nej rec; nactu rec a spocitam
1025  // sbalg a matrix se zadavaji v souboru
1026  if (1) {
1027  p = new Problem;
1028  p->read_input_file (vtuGeomFileField);
1029  p->input_data_initialize_and_check_consistency();
1030  p->convert_to_equivalent_problem ();
1031  p->print_equivalent_problem (vtkAmndFile);
1032  delete p;
1033 
1034  p = new Problem;
1035  p->read_input_file (vtkAmndFile);
1036  p->input_data_initialize_and_check_consistency();
1037  p->set_diffType (DT_ANALITICAL);
1038  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_FULL);
1039  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_FULL);
1040  delete p;
1041 
1042  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1043  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1044  // compare and delete files
1045  status += cmp_vtk_file_delete_first (vtkAmndFile);
1046  }
1047  // nactu core file a spocitam
1048  if (1) {
1049  p = new Problem;
1050  p->read_input_file (vtuGeomFileField);
1051  p->input_data_initialize_and_check_consistency();
1052  p->convert_to_equivalent_problem ();
1053  p->set_diffType (DT_ANALITICAL);
1054  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_FULL);
1055  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_FULL);
1056  delete p;
1057 
1058  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1059  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1060  }
1061  // nactu core file, vytisknu rec a spocitam
1062  if (1) {
1063  p = new Problem;
1064  p->read_input_file (vtuGeomFileField);
1065  p->input_data_initialize_and_check_consistency();
1066  p->convert_to_equivalent_problem ();
1067  p->print_equivalent_problem (vtkAmndFile);
1068  p->set_diffType (DT_ANALITICAL);
1069  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_FULL);
1070  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_FULL);
1071  delete p;
1072 
1073  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1074  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1075  // compare and delete files
1076  status += cmp_vtk_file_delete_first (vtkAmndFile);
1077  }
1078  // nactu core file, spocitam a vytisknu rec
1079  if (1) {
1080  p = new Problem;
1081  p->read_input_file (vtuGeomFileField);
1082  p->input_data_initialize_and_check_consistency();
1083  p->convert_to_equivalent_problem ();
1084  p->set_diffType (DT_ANALITICAL);
1085  p->giveFieldsOfPointOneRS(d1, e1, s1, coords1, 'p', 0, PFCM_FULL);
1086  p->giveFieldsOfPointOneRS(d2, e2, s2, coords2, 'p', 0, PFCM_FULL);
1087  p->print_equivalent_problem (vtkAmndFile);
1088  delete p;
1089 
1090  if (diff_fields (numDisp, numStrain, numStress, d1ok, d1, e1ok, e1, s1ok, s1, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1091  if (diff_fields (numDisp, numStrain, numStress, d2ok, d2, e2ok, e2, s2ok, s2, 1e-7)) { _warningg("Ellipsoid results differs"); status++; }
1092  // compare and delete files
1093  status += cmp_vtk_file_delete_first (vtkAmndFile);
1094  }
1095 
1096  return status;
1097 }
1098 
1100 int X2D_3incl (void)
1101 {
1102  int status = 0;
1103  char vtkGeomFile[] = "verified/X2D_3incl.geom.vtk";
1104  char vtkMeshFile[] = "verified/X2D_3incl.mesh.vtk";
1105  char vtkMeshFileOut[] = "verified/X2D_3incl.mesh.out";
1106  char vtkMeshFileOut0[] = "verified/X2D_3incl.mesh.out.00.vtk";
1107  char vtkMeshFileOut1[] = "verified/X2D_3incl.mesh.out.01.vtk";
1108 
1109  // compute
1110  Problem *p = new Problem;
1111  p->set_SBAM(SBAT_ORIGINAL);
1112  p->set_SBAM(SBAT_ORIGINAL);
1113  p->read_inclusions_plus_initialize_and_print (vtkGeomFile, NULL, DT_NUMERICAL);
1114  p->set_intFieldsShape(4);
1115  p->printFieldsOnMeshVTK (vtkMeshFileOut, vtkMeshFile, 'p', 0, 0, PFCM_OPTIMIZED);
1116  delete p;
1117 
1118  // compare files
1119  status += cmp_vtk_file_delete_first (vtkMeshFileOut0);
1120  status += cmp_vtk_file_delete_first (vtkMeshFileOut1);
1121  // delete files
1122  return status;
1123 }
1125 int X3D_3incl (void)
1126 {
1127  int status = 0;
1128  char vtkGeomFile[] = "verified/X3D_3incl.geom.vtk";
1129  //char vtkTopoFile2d[] = "verified/X3D_3incl.topo2d.vtk";
1130  char vtkTopoFile3d[] = "verified/X3D_3incl.topo3d.vtk";
1131  char vtkMeshFile[] = "verified/X3D_3incl.mesh.vtk";
1132  char vtkMeshFileOut[] = "verified/X3D_3incl.mesh.out";
1133  char vtkMeshFileOut0[] = "verified/X3D_3incl.mesh.out.00.vtk";
1134 
1135  // compute
1136  Problem *p = new Problem;
1137  p->set_SBAM(SBAT_ORIGINAL);
1138 
1139  p->read_inclusions_plus_initialize_and_print (vtkGeomFile, NULL, DT_ANALITICAL);
1140  p->set_intFieldsShape(4);
1141  p->printFieldsOnMeshVTK (vtkMeshFileOut, vtkMeshFile, 'p', 0, 0, PFCM_OPTIMIZED);
1142  //p->print_visualization(vtkTopoFile2d, 10, 2);
1143  p->print_visualization(vtkTopoFile3d, 10, 3);
1144  delete p;
1145 
1146  // compare and delete files
1147  status += cmp_vtk_file_delete_first (vtkMeshFileOut0);
1148  return status;
1149 }
1150 
1151 
1152 
1153 
1154 
1155 // /// 1x 3D inkluze, simulujici 2d problem
1156 // /// integrace vsech poli na regularni siti
1157 // int integral_fields (void)
1158 // {
1159 // int status = 0;
1160 // //
1161 // Problem *p = new Problem;
1162 //
1163 // // ************** DATA INPUT ******************
1164 // p->set_data_set();
1165 // p->set_dimension(3);
1166 // p->set_diffType (DT_ANALITICAL);
1167 // p->set_matrix_E_nu(1.0, 0.4);
1168 // //
1169 // p->set_numberOfRemoteStrains(1);
1170 // double rs[] = { 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0 }; p->set_RemoteStrain(0, rs);
1171 //
1172 // // loop over all inclusions
1173 // p->set_number_of_inclusions(1);
1174 // // 1. inclusion
1175 // p->set_inclusion_centroids(0,0.0, 0.0, 0.0);
1176 // p->give_inclusion(0)->set_Inclusion_shape(IS_PROLATE_SPHEROID);
1177 // p->set_inclusion_E_nu(0, 20, 0.4);
1178 // p->set_inclusion_semiaxesDimensions(0, 1000, 0.25, 0.25);
1179 // p->set_inclusion_EulerAnglesDEG(0, 0.0, 90.0, 90.0);
1180 //
1181 // p->input_data_initialize_and_check_consistency();
1182 // // ************** DATA INPUT ******************
1183 //
1184 // p->convert_to_equivalent_problem ();
1185 //
1186 // double p1[3] = { -1.0, -1.0, 0.0};
1187 // double p2[3] = { 1.0, 1.0, 0.0};
1188 // long n[3] = { 80, 80, 0};
1189 //
1190 // double *strain = new double[6]; double strainOK = -0.074459048068787;
1191 // double *stress = new double[6]; double stressOK = 0.057726312495177;
1192 //
1193 // Mesh *m = p->give_new_mesh();
1194 // m->generate_regular_mesh (p1, p2, n);
1195 // m->compute_element_fields ('p',true, true, true, 0, 1, PFCM_OPTIMIZED);
1196 // m->integrate_pertFields(0, &strain, &stress, 0, 1);
1197 //
1198 // if (compare_fields (1, strain+3, &strainOK, 1e-7)) { _warningg("integral_fields results differs"); status++; }
1199 // if (compare_fields (1, stress+3, &stressOK, 1e-7)) { _warningg("integral_fields results differs"); status++; }
1200 //
1201 // delete [] strain;
1202 // delete [] stress;
1203 // delete p;
1204 // return status;
1205 // }
1206 
1207 
1208 /*end of file*/
x3D_1I_Ellipsoid
int x3D_1I_Ellipsoid(void)
Single inclusion of type ellipsoid.
Definition: tests_verified.cpp:318
plstrain_2I_3d
int plstrain_2I_3d(void)
Two inclusions of type ellipse by 3d.
Definition: tests_verified.cpp:532
mumech::Problem::input_data_initialize_and_check_consistency
void input_data_initialize_and_check_consistency(void)
Initializes and checks consistency of all input data. This function has to be called after data input...
Definition: problem.cpp:127
grid_2d_3x3_3d_F
int grid_2d_3x3_3d_F(void)
Nine 3d inclusions in 2d grid 3x3.
Definition: tests_verified.cpp:405
_warningg
#define _warningg(_1)
Definition: gelib.h:163
mumech::Problems::set_dimension
void set_dimension(int d)
Definition: problem.h:109
types.h
file of various types and symbolic constant definitions
x2D_1I_Ellipse_homog
int x2D_1I_Ellipse_homog(void)
Single inclusion of type ellipse.
Definition: tests_verified.cpp:604
mumech::Problem::give_nLC
virtual int give_nLC(void) const
Give number of load cases, e.i. number of load cases, i.e. number of Remote_strain fields...
Definition: problem.h:334
mumech::Problem::set_inclusion_semiaxesDimensions
void set_inclusion_semiaxesDimensions(long id, double a1, double a2, double a3=0.0)
Sets the id -th inclusion sorted semiaxes dimensions [a1,a2,a3]. The a3 dimension is ignored in the c...
Definition: problem.h:440
mumech::Problem::read_inclusions_plus_initialize_and_print
Problem * read_inclusions_plus_initialize_and_print(const char *inclusion_file, const char *equiv_file, DiffTypes dt)
Definition: problem.cpp:2059
mumech::Problem
Problem description.
Definition: problem.h:154
mumech::SBAT_ORIGINAL
Original Honza Novak&#39;s balancing.
Definition: types.h:90
Ellipse
int Ellipse(void)
Definition: tests_verified.cpp:784
librw.h
Input / output function.
mumech::Inclusion::give_StiffnessMatrixFull
void give_StiffnessMatrixFull(double *answer) const
Copy stiffness tensor into full matrix answer.
Definition: inclusion.cpp:799
mumech::Problem::set_inclusion_E_nu
void set_inclusion_E_nu(long id, double E, double nu)
Sets the id -th inclusion material properties, Young&#39;s modulus E and Poissons ratio nu...
Definition: problem.h:438
mumech::Problem::print_equivalent_problem
void print_equivalent_problem(const char *filename)
Prints the equivalent problem record into the VTK file.
Definition: problem.cpp:2122
mumech::Homogenization
Class of function for homogenization of stress fields.
Definition: homogenization.h:46
mesh.h
Class Mesh.
mumech::Problem::give_inclusion
Inclusion * give_inclusion(long i) const
debug, for one function in tools
Definition: problem.h:402
mumech::Problem::set_inclusion_EulerAnglesDEG
void set_inclusion_EulerAnglesDEG(long id, double e1, double e2=0.0, double e3=0.0)
Sets the id -th inclusion Euler angles [e1,e2,e3] given in degrees. The e2 and e3 angles are ignored ...
Definition: problem.h:442
melement.h
Class mElement, mesh element.
tests_verified
void tests_verified(void)
Definition: tests_verified.cpp:49
grid_2d_3x3_3d
int grid_2d_3x3_3d(void)
Nine 3d inclusions in 2d grid 3x3.
Definition: tests_verified.cpp:351
mumech::Problem::set_RemoteStrain
void set_RemoteStrain(int id, const double *rs)
Sets the id -th remote strain, rs is row-by-row matrix of dimensions 2x2 or 3x3 for 2d or 3d problem...
Definition: problem.h:431
compare_fields
int compare_fields(int n, const double *f1, const double *f2, double eps)
Compare two fields.
Definition: tests_verified.cpp:85
mumech::Problem::set_inclusion_centroids
void set_inclusion_centroids(long id, double x, double y, double z=0.0)
Sets the id -th inclusion centroids [x,y,z]. The z coordinate is ignored in the case of 2 dimensions...
Definition: problem.h:436
diff_fields
int diff_fields(int numDisp, int numStrain, int numStress, const double *d1, const double *d2, const double *e1, const double *e2, const double *s1, const double *s2, double eps)
Definition: tests_verified.cpp:97
mumech::Homogenization::set_boundingBox
void set_boundingBox(double x1, double y1, double x2, double y2)
Definition: homogenization.cpp:62
x2D_1I_Ellipse_direct_API
int x2D_1I_Ellipse_direct_API(void)
Single inclusion of type ellipse.
Definition: tests_verified.cpp:165
mumech::Problem::set_SBA_optimized
void set_SBA_optimized(bool val)
Set type of self-balancing algorithm.
Definition: problem.h:368
x3D_1I_Ellipsoid_homog
int x3D_1I_Ellipsoid_homog(void)
Single inclusion of type ellipse.
Definition: tests_verified.cpp:662
cmp_txt_file_delete_first
int cmp_txt_file_delete_first(const char *file1, const char *file2)
Definition: tests_verified.cpp:135
mumech::Problem::matrix_giveFullStiffMatrix
void matrix_giveFullStiffMatrix(double *m) const
Definition: problem.h:330
mumech::Problem::set_number_of_inclusions
void set_number_of_inclusions(long n)
Sets number of inclusions.
Definition: problem.cpp:136
mumech::Problem::set_SBAM
void set_SBAM(SBAtype val)
Definition: problem.h:372
cmp_vtk_file_delete_first
int cmp_vtk_file_delete_first(const char *file1, const char *file2)
Definition: tests_verified.cpp:110
comparison.h
Class Comparison.
X3D_3incl
int X3D_3incl(void)
Definition: tests_verified.cpp:1125
mumech::Problem::set_data_set
void set_data_set(void)
Function switches bool flag data_set on true. That indicates start of the process of data input...
Definition: problem.h:423
mumech::Problem::set_matrix_E_nu
void set_matrix_E_nu(double E, double nu)
Sets problem dimension.
Definition: problem.h:427
mumech::Problem::set_DataDimDiff
void set_DataDimDiff(int dim, DiffTypes dt)
Function agregates set_data_set, set_dimension and set_diffType functions.
Definition: problem.h:457
mumech::Homogenization::giveTotalVolumeFractionOfInclusions
double giveTotalVolumeFractionOfInclusions(void) const
Definition: homogenization.cpp:172
homogenization.h
Class Homogenization.
gelibspace::FP_cmp_files
bool FP_cmp_files(const char *file1, const char *file2)
Return false if no differences.
Definition: librw.cpp:230
vtk.h
I/O VTK functions.
mumech::Problem::set_diffType
void set_diffType(DiffTypes val)
Sets type of differentiation, analytical by default.
Definition: problem.h:455
mumech::Problem::give_new_homogenization
Homogenization * give_new_homogenization(HomogenizationType ht)
Definition: problem.cpp:2084
Four_3d_inclusions
int Four_3d_inclusions(void)
Four 3d inclusion.
Definition: tests_verified.cpp:727
Ellipsoid
int Ellipsoid(void)
Definition: tests_verified.cpp:866
X2D_3incl
int X2D_3incl(void)
Definition: tests_verified.cpp:1100
gelibspace::file_delete_if_exist
bool file_delete_if_exist(const char *name)
Definition: gelib.h:427
mumech::Problem::set_inclusion_all
void set_inclusion_all(long id, double x, double y, double e, double n, double a1, double a2, double e1)
Agregates ...
Definition: problem.h:464
mumech::Problem::printFieldsOnMeshGrid
void printFieldsOnMeshGrid(const char *mesh_file_out, const double *p1, const double *p2, const long *n, char ptFlag, int rs, int nrs, PFCmode pfcMode=PFCM_OPTIMIZED) const
Function computes all fields (displacements, strains and stresses) in nodes a regular orthogonal mesh...
Definition: problem.cpp:2041
mumech::Homogenization::giveHomogenizedStiffnessMatrix
virtual void giveHomogenizedStiffnessMatrix(double *answer)=0
Function returning the homogenized stiffness matrix according the defined method termitovo pridat par...
mumech::Problem::print_visualization
void print_visualization(const char *filename, int n, int dim=0, bool refined=false)
Triangulates inclusions surfaces and prints filename VTK file.
Definition: problem.cpp:2237
mumech::Problem::read_input_file
void read_input_file(const char *filename)
Reads filename VTK input file containing inclusion geometry and topology.
Definition: problem.cpp:357
mumech::cmp_vtk_file
bool cmp_vtk_file(const char *file1, const char *file2, int verbose)
Compare two VTK files. Return false if same. [ return hodnota byly kdysi opacna ].
Definition: vtk.cpp:268
mumech::Problem::printFieldsOnMeshVTK
void printFieldsOnMeshVTK(const char *mesh_file_out, const char *mesh_file, char ptFlag, int rs, int nrs, PFCmode pfcMode=PFCM_OPTIMIZED) const
Function computes all fields (displacements, strains and stresses) in nodes a mesh given via VTK file...
Definition: problem.cpp:2028
gelibspace::file_delete_if_exist_exit
void file_delete_if_exist_exit(const char *name)
Definition: gelib.h:437
mumech::Homogenization::giveClassName
virtual const char * giveClassName()
Function returning class name.
Definition: homogenization.h:66
x3D_1I_Ellipsoid_direct_API
int x3D_1I_Ellipsoid_direct_API(void)
Single inclusion of type ellipsoid.
Definition: tests_verified.cpp:239
mumech::Problem::set_intFieldsShape
void set_intFieldsShape(int val)
Definition: problem.h:371
plstrain_2I_2d
int plstrain_2I_2d(void)
Two inclusions of type ellipse.
Definition: tests_verified.cpp:459
mumech::Problem::set_UnitRemoteStrains
void set_UnitRemoteStrains(void)
Function sets a set of unit remote strains.
Definition: problem.h:462
problem.h
Class Problem.
mumech::Problem::set_numberOfRemoteStrains
void set_numberOfRemoteStrains(int n)
Sets number of remote strains.
Definition: problem.h:429
mumech::Problem::giveFieldsOfPointOneRS
long giveFieldsOfPointOneRS(double *displc, double *strain, double *stress, const double *coords, char ptFlag, int rs, PFCmode pfcMode=PFCM_OPTIMIZED, long reqIncl=-3, STRNotation tn=STRN_THEORETICAL_ROW) const
Function gives results for only one remote strain. See function giveFieldsOfPoint().
Definition: problem.cpp:1699
mumech::e
void e(int i)
Definition: matrix_vector.cpp:45
mumech::Problem::convert_to_equivalent_problem
void convert_to_equivalent_problem(void)
Converts the given heterogeneous problem to the equivalent problem.
Definition: problem.cpp:781
This reference manual is part of the muMECH documentation. Generated on Wed Feb 24 2016 22:50:08 for muMECH by doxygen 1.8.11.