CArl
Code Arlequin / C++ implementation
libmesh_assemble_lin_anisotropic.cpp File Reference

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 3 of file libmesh_assemble_lin_anisotropic.cpp.

3  {
4 
5  // --- Initialize libMesh
6  libMesh::LibMeshInit init(argc, argv);
7 
8  // Do performance log?
9  const bool MASTER_bPerfLog_carl_libmesh = true;
10  libMesh::PerfLog perf_log("Main program", MASTER_bPerfLog_carl_libmesh);
11 
12  // libMesh's C++ / MPI communicator wrapper
13  libMesh::Parallel::Communicator& WorldComm = init.comm();
14 
15  // Number of processors and processor rank.
16  int rank = WorldComm.rank();
17  int nodes = WorldComm.size();
18 
19  // --- Set up inputs
20 
21  // Command line parser
22  GetPot command_line(argc, argv);
23 
24  // File parser
25  GetPot field_parser;
26 
27  // If there is an input file, parse it to get the parameters. Else, parse the command line
28  std::string input_filename;
29  if (command_line.search(2, "--inputfile", "-i")) {
30  input_filename = command_line.next(input_filename);
31  field_parser.parse_input_file(input_filename, "#", "\n", " \t\n");
32  } else {
33  field_parser = command_line;
34  }
35 
36  libmesh_assemble_input_params input_params;
37  get_input_params(field_parser, input_params);
38 
39  // Check libMesh installation dimension
40  const unsigned int dim = 3;
41 
42  libmesh_example_requires(dim == LIBMESH_DIM, "3D support");
43 
44  // --- Declare the three meshes to be intersected
45 
46  // - Parallelized meshes: A, B, mediator and weight
47  perf_log.push("Meshes - Parallel","Read files:");
48  libMesh::Mesh system_mesh(WorldComm, dim);
49  system_mesh.read(input_params.mesh_file);
50  system_mesh.prepare_for_use();
51 
52  libMesh::Mesh mesh_weight(WorldComm, dim);
53  mesh_weight.allow_renumbering(false);
54  mesh_weight.read(input_params.mesh_weight_file);
55  mesh_weight.prepare_for_use();
56 
57  perf_log.pop("Meshes - Parallel","Read files:");
58 
59 // // --- Generate the equation systems
60 // perf_log.push("Initialization","System initialization:");
61 // carl::coupled_system CoupledTest(WorldComm,input_params.solver_type);
62 
63 // libMesh::EquationSystems& equation_systems_inter =
64 // CoupledTest.add_inter_EquationSystem("InterSys", mesh_inter);
65 
66 // // Add the weight function mesh
67 // CoupledTest.add_alpha_mask("MicroSys",mesh_weight);
68 // CoupledTest.set_alpha_mask_parameters("MicroSys",domain_Idx_BIG,domain_Idx_micro[0],domain_Idx_coupling[0]);
69 
70 // perf_log.pop("Initialization","System initialization:");
71 
72 // // - Build the MACRO system
73 
74 // perf_log.push("Macro system","System initialization:");
75 // libMesh::EquationSystems& equation_systems_BIG =
76 // CoupledTest.set_BIG_EquationSystem("BigSys", system_mesh);
77 
78 // // [MACRO] Set up the physical properties
79 // libMesh::LinearImplicitSystem& elasticity_system_BIG
80 // = add_elasticity(equation_systems_BIG);
81 
82 // // [MACRO] Defining the boundaries with Dirichlet conditions
83 // set_displaced_border_translation(elasticity_system_BIG, z_upper_hole,z_upper_hole_id);
84 // set_displaced_border_translation(elasticity_system_BIG, z_lower_hole,z_lower_hole_id);
85 // set_clamped_border(elasticity_system_BIG, fixed_bound_id);
86 
87 // // [MACRO] Build stress system
88 // libMesh::ExplicitSystem& stress_system_BIG
89 // = add_stress(equation_systems_BIG);
90 
91 // equation_systems_BIG.init();
92 
93 // perf_log.pop("Macro system","System initialization:");
94 
95 // // - Build the micro system
96 
97 // perf_log.push("Micro system","System initialization:");
98 // libMesh::EquationSystems& equation_systems_micro =
99 // CoupledTest.add_micro_EquationSystem<libMesh::PetscMatrix<libMesh::Number> >("MicroSys", mesh_micro);
100 
101 // // [MICRO] Set up the physical properties
102 // libMesh::LinearImplicitSystem& elasticity_system_micro
103 // = add_elasticity(equation_systems_micro);
104 
105 // // [MICRO] Build stress system
106 // libMesh::ExplicitSystem& stress_system_micro
107 // = add_stress(equation_systems_micro);
108 
109 // equation_systems_micro.init();
110 // perf_log.pop("Micro system","System initialization:");
111 
112 // // - Build the RESTRICTED BIG system
113 // perf_log.push("RESTRICTED macro system","System initialization:");
114 // libMesh::EquationSystems& equation_systems_R_BIG =
115 // CoupledTest.set_Restricted_BIG_EquationSystem("BigSys", mesh_R_BIG);
116 
117 // // [R. MACRO] Set up the physical properties
118 // libMesh::ExplicitSystem& elasticity_system_R_BIG
119 // = add_explicit_elasticity(equation_systems_R_BIG);
120 
121 // carl::reduced_system_init(elasticity_system_R_BIG);
122 
123 // perf_log.pop("RESTRICTED macro system","System initialization:");
124 
125 // // - Build the RESTRICTED micro system
126 
127 // perf_log.push("RESTRICTED micro system","System initialization:");
128 // libMesh::EquationSystems& equation_systems_R_micro =
129 // CoupledTest.add_Restricted_micro_EquationSystem("MicroSys", mesh_R_micro);
130 
131 // // [R. MICRO] Set up the physical properties
132 // libMesh::ExplicitSystem& elasticity_system_R_micro
133 // = add_explicit_elasticity(equation_systems_R_micro);
134 
135 // carl::reduced_system_init(elasticity_system_R_micro);
136 // perf_log.pop("RESTRICTED micro system","System initialization:");
137 
138 // // - Build the mediator system
139 
140 // perf_log.push("Mediator system","System initialization:");
141 // libMesh::EquationSystems& equation_systems_mediator =
142 // CoupledTest.add_mediator_EquationSystem("MediatorSys", mesh_mediator);
143 
144 // libMesh::LinearImplicitSystem& elasticity_system_mediator
145 // = add_elasticity(equation_systems_mediator);
146 
147 // equation_systems_mediator.init();
148 
149 // WorldComm.barrier();
150 
151 // perf_log.pop("Mediator system","System initialization:");
152 
153 // // - Build the dummy inter system
154 
155 // perf_log.push("Intersection system","System initialization:");
156 // libMesh::ExplicitSystem& elasticity_system_inter
157 // = add_explicit_elasticity(equation_systems_inter);
158 
159 // carl::reduced_system_init(elasticity_system_inter);
160 
161 // WorldComm.barrier();
162 
163 // perf_log.pop("Intersection system","System initialization:");
164 
165 // perf_log.push("Physical properties","System initialization:");
166 
167 // // --- Set the physical parameters
168 // double BIG_E = 0;
169 // double BIG_Mu = 0;
170 
171 // double coupling_const = -1;
172 
173 // // Anisotropic properties for the micro system
174 // carl::anisotropic_elasticity_tensor_cubic_sym anisotropy_data(equation_systems_micro,input_params.physical_params_file,BIG_E,BIG_Mu);
175 
176 // // Homogeneous properties for the macro system
177 // set_constant_physical_properties(equation_systems_BIG,BIG_E,BIG_Mu);
178 
179 // perf_log.pop("Physical properties","System initialization:");
180 
181 // // --- Set the coupling matrices
182 // perf_log.push("Set coupling matrices");
183 // coupling_const = BIG_E;
184 
185 // // Set parameters
186 // CoupledTest.set_coupling_parameters("MicroSys",coupling_const,input_params.mean_distance);
187 
188 // // Use H1 norm
189 // CoupledTest.use_H1_coupling("MicroSys");
190 
191 // // Assemble!
192 // CoupledTest.assemble_coupling_elasticity_3D_parallel("BigSys","MicroSys",
193 // "InterSys","MediatorSys",
194 // mesh_R_BIG, mesh_R_micro,
195 // local_intersection_pairs_map,
196 // local_intersection_restricted_pairs_map,
197 // local_intersection_meshI_to_inter_map,
198 // inter_mediator_BIG,
199 // inter_mediator_micro);
200 
201 // // Print some info
202 // std::cout << std::endl;
203 // std::cout << "| ---> Constants " << std::endl;
204 // std::cout << "| Macro :" << std::endl;
205 // std::cout << "| E : " << BIG_E << std::endl;
206 // std::cout << "| Mu (lamba_2) : " << BIG_Mu << std::endl;
207 // std::cout << "| lambda_1 : " << eval_lambda_1(BIG_E,BIG_Mu) << std::endl;
208 // std::cout << "| Coupling :" << std::endl;
209 // std::cout << "| kappa : " << coupling_const << std::endl;
210 // std::cout << "| e : " << input_params.mean_distance << std::endl;
211 
212 // switch(input_params.solver_type)
213 // {
214 // case carl::LATIN_MODIFIED_STIFFNESS:
215 // case carl::LATIN_ORIGINAL_STIFFNESS:
216 // {
217 // std::cout << "| LATIN :" << std::endl;
218 // std::cout << "| k_dA, k_dB : " << input_params.k_dA << " " << input_params.k_dB << std::endl;
219 // std::cout << "| k_cA, k_cB : " << input_params.k_cA << " " << input_params.k_cB << std::endl;
220 // break;
221 // }
222 // case carl::CG:
223 // {
224 // break;
225 // }
226 // }
227 // perf_log.pop("Set coupling matrices");
228 
229 // std::cout << "| restart file : " << input_params.coupled_restart_file_base << "*" << std::endl;
230 
231 // std::cout << std::endl << "| --> Testing the solver " << std::endl << std::endl;
232 // perf_log.push("Set up","Coupled Solver:");
233 
234 // // --- Set restart files
235 // if(input_params.b_PrintRestartFiles || input_params.b_UseRestartFiles)
236 // {
237 // CoupledTest.set_restart( input_params.b_UseRestartFiles,
238 // input_params.b_PrintRestartFiles,
239 // input_params.coupled_restart_file_base);
240 // }
241 // std::cout << std::endl << "| --> Setting the solver " << std::endl << std::endl;
242 
243 // // --- Set up the coupled solver!
244 // switch(input_params.solver_type)
245 // {
246 // case carl::LATIN_MODIFIED_STIFFNESS:
247 // case carl::LATIN_ORIGINAL_STIFFNESS:
248 // {
249 // CoupledTest.set_LATIN_solver( "MicroSys","Elasticity",
250 // assemble_elasticity_with_weight,
251 // assemble_elasticity_anisotropic_with_weight,
252 // anisotropy_data,
253 // input_params.k_dA, input_params.k_dB, input_params.k_cA, input_params.k_cB,
254 // input_params.LATIN_eps, input_params.LATIN_conv_max, input_params.LATIN_relax);
255 // break;
256 // }
257 // case carl::CG:
258 // {
259 // CoupledTest.use_null_space_micro("MicroSys",true);
260 // CoupledTest.set_cg_preconditioner_type(input_params.CG_precond_type);
261 // CoupledTest.set_CG_solver( "MicroSys","Elasticity",
262 // assemble_elasticity_with_weight,
263 // assemble_elasticity_anisotropic_with_weight,
264 // anisotropy_data,
265 // input_params.CG_coupled_conv_abs,input_params.CG_coupled_conv_rel,
266 // input_params.CG_coupled_conv_max,input_params.CG_coupled_div,
267 // input_params.CG_coupled_conv_corr);
268 // break;
269 // }
270 // }
271 // perf_log.pop("Set up","Coupled Solver:");
272 
273 
274 // // --- Solve !
275 // perf_log.push("Solve","Coupled Solver:");
276 // CoupledTest.solve("MicroSys","Elasticity",input_params.coupled_convergence_output);
277 // perf_log.pop("Solve","Coupled Solver:");
278 
279 // // --- Calculate stress
280 // perf_log.push("Compute stress - micro","Output:");
281 // compute_stresses(equation_systems_micro);
282 // perf_log.pop("Compute stress - micro","Output:");
283 
284 // perf_log.push("Compute stress - macro","Output:");
285 // compute_stresses(equation_systems_BIG);
286 // perf_log.pop("Compute stress - macro","Output:");
287 
288 // // --- Export solution
289 // #ifdef LIBMESH_HAVE_EXODUS_API
290 // if(input_params.b_PrintOutput)
291 // {
292 // perf_log.push("Save output","Output:");
293 // libMesh::ExodusII_IO exo_io_micro(mesh_micro, /*single_precision=*/true);
294 
295 // std::set<std::string> system_names_micro;
296 // system_names_micro.insert("Elasticity");
297 // exo_io_micro.write_equation_systems(input_params.output_file_micro,equation_systems_micro,&system_names_micro);
298 
299 // exo_io_micro.write_element_data(equation_systems_micro);
300 
301 // libMesh::ExodusII_IO exo_io_BIG(system_mesh, /*single_precision=*/true);
302 
303 // std::set<std::string> system_names_BIG;
304 // system_names_BIG.insert("Elasticity");
305 // exo_io_BIG.write_equation_systems(input_params.output_file_BIG,equation_systems_BIG,&system_names_BIG);
306 
307 // exo_io_BIG.write_element_data(equation_systems_BIG);
308 // perf_log.pop("Save output","Output:");
309 // }
310 // #endif
311 
312 // if(input_params.b_ExportScalingData)
313 // {
314 // CoupledTest.print_perf_log(input_params.scaling_data_file);
315 // }
316 
317  return 0;
318 }
std::string mesh_weight_file
Path to the mesh containing the weight region indices.
void get_input_params(GetPot &field_parser, feti_iterate_params &input_params)
Parser function for the coupled solver test programs.
std::string mesh_file
Path to the system mesh.