CArl
Code Arlequin / C++ implementation
assemble_functions_coupling.cpp
Go to the documentation of this file.
1 #include "assemble_coupling.h"
2 
3 namespace carl
4 {
5 
6  void assemble_coupling_matrices::prepare_coupling_preallocation(
7  libMesh::PetscMatrix<libMesh::Number>& coupling_matrix,
8  libMesh_fe_addresses_3& row_addresses,
9  libMesh_fe_addresses_3& col_addresses,
10  const std::unordered_multimap<int,int>& inter_table
11  )
12 {
13  const libMesh::Parallel::Communicator& WorldComm = row_addresses.mesh.comm();
14  int rank = WorldComm.rank();
15  int nodes = WorldComm.size();
16 
17  // Local and global dimensions setup
18  const unsigned int row_M = row_addresses.dof_map.n_dofs();
19  const int col_N = col_addresses.dof_map.n_dofs();
20  const unsigned int row_M_local = row_addresses.dof_map.n_local_dofs();
21  const int col_N_local = col_addresses.dof_map.n_local_dofs();
22 
23  // Vectors that will be used to calculate the preallocation. There are two
24  // "versions": a local one, starting at zero and which will contain this
25  // processor's preallocation after syncing, and a global one, which will be
26  // send to other processors
27  std::vector<unsigned int> local_n_nz(row_M_local,0);
28  std::vector<unsigned int> local_n_oz(row_M_local,0);
29 
30  // Local limits of the DoF's
31  std::vector<int> begin_row_DoF(nodes);
32  std::vector<int> end_row_DoF(nodes);
33  for(int iii = 0; iii < nodes; ++iii)
34  {
35  begin_row_DoF[iii] = row_addresses.dof_map.first_dof(iii);
36  end_row_DoF[iii] = row_addresses.dof_map.end_dof(iii);
37  }
38 
39  std::vector<int> begin_col_DoF(nodes);
40  std::vector<int> end_col_DoF(nodes);
41  for(int iii = 0; iii < nodes; ++iii)
42  {
43  begin_col_DoF[iii] = col_addresses.dof_map.first_dof(iii);
44  end_col_DoF[iii] = col_addresses.dof_map.end_dof(iii);
45  }
46 
47  int local_begin_row_DoF = begin_row_DoF[rank];
48  int local_end_row_DoF = end_row_DoF[rank];
49 
50  // Iterate over the row's mesh cells
51 
52  int row_mesh_elem_idx = 0;
53  int dummy_col_dof;
54 
55  std::vector<std::vector<int> > to_add_per_proc(row_M, std::vector<int>(nodes,0));
56 
57  std::vector<int> row_dof_indices_procs;
58  libMesh::MeshBase::const_element_iterator it_row_elem =
59  row_addresses.mesh.active_local_elements_begin();
60  const libMesh::MeshBase::const_element_iterator it_row_elem_end =
61  row_addresses.mesh.active_local_elements_end();
62 
63  // What the code has to do:
64  /*
65  * - For each element associated to the rows, get the intersecting
66  * elements and find out which processors own each of its DoFs. Use
67  * this to populate the "to_add_per_proc" table.
68  *
69  * - Once it is filled, sync the "to_add_per_proc" between the
70  * processors.
71  *
72  * - On each processor, iterate over the local row DoFs and add the
73  * info from "to_add_per_proc" to either "local_n_nz" or "local_n_oz",
74  * depending on the processor.
75  *
76  * - ... pray that it works?
77  */
78 
79  std::unordered_set<int> col_DoFs_to_allocate;
80  col_DoFs_to_allocate.reserve(col_addresses.mesh.n_nodes());
81 
82  for( ; it_row_elem != it_row_elem_end; ++it_row_elem )
83  {
84  // For each local mediator mesh elements ...
85  const libMesh::Elem* elem_mediator = *it_row_elem;
86 
87  // ... extract the element's index ...
88  row_mesh_elem_idx = elem_mediator->id();
89 
90  // ... set its DoFs ...
91  row_addresses.set_DoFs(row_mesh_elem_idx);
92 
93  // ... clear the col_DoFs set
94  col_DoFs_to_allocate.clear();
95 
96  // ... get the DoFs from the col mesh ...
97  auto range_elem_col_idx = inter_table.equal_range(row_mesh_elem_idx);
98 
99  for(auto it = range_elem_col_idx.first; it != range_elem_col_idx.second; ++it)
100  {
101  col_addresses.set_DoFs(it->second);
102 
103  for(unsigned int iii = 0; iii < col_addresses.n_dofs; ++iii)
104  {
105  col_DoFs_to_allocate.insert(col_addresses.dof_indices[iii]);
106  }
107  }
108 
109  // ... and finally increase the preallocation
110  for(auto it = col_DoFs_to_allocate.begin(); it != col_DoFs_to_allocate.end(); ++it)
111  {
112  dummy_col_dof = *it;
113 
114  for(int nnn = 0; nnn < nodes; ++nnn)
115  {
116  if( dummy_col_dof > begin_col_DoF[nnn] &&
117  dummy_col_dof < end_col_DoF[nnn])
118  {
119  for(unsigned int jjj = 0; jjj < row_addresses.dof_indices.size(); ++jjj)
120  {
121  ++to_add_per_proc[row_addresses.dof_indices[jjj]][nnn];
122  }
123  break;
124  }
125  }
126  }
127  }
128 
129  // Guarantee that everyone is in the same page
130  WorldComm.barrier();
131 
132  // Do a (lot) of syncs
133  for(int nnn = 0; nnn < nodes; ++nnn)
134  {
135  for(int iii = begin_row_DoF[nnn]; iii < end_row_DoF[nnn]; ++iii)
136  {
137  MPI_reduce_vector(to_add_per_proc[iii],nnn,WorldComm);
138  }
139  }
140 
141  int dummy_total_nz = 0;
142  int dummy_total_oz = 0;
143 
144  // Now, calculate the preallocation
145  for(int iii = local_begin_row_DoF; iii < local_end_row_DoF; ++iii)
146  {
147  for(int nnn = 0; nnn < rank; ++nnn)
148  {
149  local_n_oz[iii - local_begin_row_DoF] += to_add_per_proc[iii][nnn];
150  dummy_total_oz += to_add_per_proc[iii][nnn];
151  }
152 
153  local_n_nz[iii - local_begin_row_DoF] += to_add_per_proc[iii][rank];
154  dummy_total_nz += to_add_per_proc[iii][rank];
155 
156  for(int nnn = rank + 1; nnn < nodes; ++nnn)
157  {
158  local_n_oz[iii - local_begin_row_DoF] += to_add_per_proc[iii][nnn];
159  dummy_total_oz += to_add_per_proc[iii][nnn];
160  }
161  }
162 
163  dummy_total_nz = 0;
164  dummy_total_oz = 0;
165 
166  for(int iii = local_begin_row_DoF; iii < local_end_row_DoF; ++iii)
167  {
168  local_n_nz[iii - local_begin_row_DoF] = 1.1*local_n_nz[iii - local_begin_row_DoF];
169  local_n_nz[iii - local_begin_row_DoF] = std::max(static_cast<unsigned int>(30),local_n_nz[iii - local_begin_row_DoF]);
170  local_n_nz[iii - local_begin_row_DoF] = std::min(static_cast<unsigned int>(col_N_local),local_n_nz[iii - local_begin_row_DoF]);
171 
172  local_n_oz[iii - local_begin_row_DoF] = 1.1*local_n_oz[iii - local_begin_row_DoF];
173  local_n_oz[iii - local_begin_row_DoF] = std::max(static_cast<unsigned int>(10),local_n_oz[iii - local_begin_row_DoF]);
174  local_n_oz[iii - local_begin_row_DoF] = std::min(static_cast<unsigned int>(col_N - col_N_local),local_n_oz[iii - local_begin_row_DoF]);
175 
176  dummy_total_nz += local_n_nz[iii - local_begin_row_DoF];
177  dummy_total_oz += local_n_oz[iii - local_begin_row_DoF];
178  }
179 
180  coupling_matrix.init(row_M, col_N, row_M_local, col_N_local, local_n_nz, local_n_oz);
181 }
182 
183  void assemble_coupling_matrices::assemble_coupling_elasticity_3D_parallel(
184  const std::string BIG_name,
185  const std::string micro_name,
186  const std::string inter_name,
187  const std::string mediator_name,
188 
189  const libMesh::MeshBase& mesh_R_BIG,
190  const libMesh::MeshBase& mesh_R_micro,
191 
192  const std::unordered_map<int,std::pair<int,int> >&
193  full_intersection_pairs_map,
194  const std::unordered_map<int,std::pair<int,int> >&
195  full_intersection_restricted_pairs_map,
196  const std::unordered_map<int,int>&
197  local_intersection_meshI_to_inter_map,
198  const std::unordered_multimap<int,int>& inter_table_mediator_BIG,
199  const std::unordered_multimap<int,int>& inter_table_mediator_micro,
200 
201  const std::string BIG_type,
202  const std::string micro_type,
203  bool bSameElemsType)
204  {
205  // TODO : make it possible to invert the algorithm's systems!
206  // Addresses to the eq. systems
207  libMesh::EquationSystems& mediator_eq_system =
208  *m_mediator_EquationSystemMap[mediator_name];
209  libMesh::EquationSystems& BIG_eq_system =
210  *m_BIG_EquationSystem.second;
211  libMesh::EquationSystems& micro_eq_system =
212  *m_micro_EquationSystemMap[micro_name];
213  libMesh::EquationSystems& inter_eq_system =
214  *m_inter_EquationSystemMap[inter_name];
215 
216  libMesh::EquationSystems& R_BIG_eq_system =
217  *m_R_BIG_EquationSystem.second;
218  libMesh::EquationSystems& R_micro_eq_system =
219  *m_R_micro_EquationSystemMap[micro_name];
220 
221  // First, test if all the systems have an elasticity model and variable set
222  homemade_assert_msg(micro_eq_system.has_system(micro_type),
223  " Micro equation systems is missing a system type!");
224  homemade_assert_msg(BIG_eq_system.has_system(BIG_type),
225  " Macro equation systems is missing a system type!");
226  homemade_assert_msg(R_micro_eq_system.has_system("Elasticity"),
227  " Restricted micro equation systems missing \"Elasticity\" system!");
228  homemade_assert_msg(R_BIG_eq_system.has_system("Elasticity"),
229  " Restricted macro equation systems missing \"Elasticity\" system!");
230  homemade_assert_msg(inter_eq_system.has_system("Elasticity"),
231  " Intersection equation systems missing \"Elasticity\" system!");
232  homemade_assert_msg(mediator_eq_system.has_system("Elasticity"),
233  " Mediatored equation systems missing \"Elasticity\" system!");
234 
235  // Systems and vars
236  libMesh::System& volume_mediator_system =
237  libMesh::cast_ref<libMesh::System&>(mediator_eq_system.get_system("Elasticity"));
238 
239  libMesh::System& volume_BIG_system =
240  libMesh::cast_ref<libMesh::System&>(BIG_eq_system.get_system<libMesh::ExplicitSystem>(BIG_type));
241 
242  libMesh::System& volume_micro_system =
243  libMesh::cast_ref<libMesh::System&>(micro_eq_system.get_system<libMesh::ExplicitSystem>(micro_type));
244 
245  libMesh::System& volume_inter_system =
246  inter_eq_system.get_system<libMesh::ExplicitSystem>("Elasticity");
247 
248  libMesh::System& volume_R_BIG_system =
249  R_BIG_eq_system.get_system<libMesh::ExplicitSystem>("Elasticity");
250 
251  libMesh::System& volume_R_micro_system =
252  R_micro_eq_system.get_system<libMesh::ExplicitSystem>("Elasticity");
253 
254  libMesh_fe_addresses_3 mediator_addresses(volume_mediator_system);
255  libMesh_fe_addresses_3 BIG_addresses(volume_BIG_system);
256  libMesh_fe_addresses_3 micro_addresses(volume_micro_system);
257  libMesh_fe_addresses_3 inter_addresses(volume_inter_system);
258 
259  libMesh_fe_addresses_3 R_BIG_addresses(volume_R_BIG_system);
260  libMesh_fe_addresses_3 R_micro_addresses(volume_R_micro_system);
261 
262  // Vector that will keep the quadrature points
263  const std::vector<libMesh::Point>& quad_points_inter =
264  inter_addresses.fe_unique_ptr->get_xyz();
265  std::vector<libMesh::Point> quad_points_reference;
266 
267  // Jacobians
268  const std::vector<libMesh::Real>& JxW =
269  inter_addresses.fe_unique_ptr->get_JxW();
270 
271  // Shape functions
272  const std::vector<std::vector<libMesh::Real> >& phi_inter =
273  inter_addresses.fe_unique_ptr->get_phi();
274  std::vector<std::vector<libMesh::Real> > corrected_phi_R_BIG;
275  std::vector<std::vector<libMesh::Real> > corrected_phi_R_micro;
276  std::vector<std::vector<libMesh::Real> > corrected_phi_mediator;
277 
278  // Shape functions gradients
279  const std::vector<std::vector<libMesh::RealGradient> >& dphi_inter =
280  inter_addresses.fe_unique_ptr->get_dphi();
281  std::vector<std::vector<libMesh::RealGradient> > corrected_dphi_R_BIG;
282  std::vector<std::vector<libMesh::RealGradient> > corrected_dphi_R_micro;
283  std::vector<std::vector<libMesh::RealGradient> > corrected_dphi_mediator;
284 
285  unsigned int n_quadrature_pts = 0;
286 
287  // Addresses to the matrices
288  libMesh::PetscMatrix<libMesh::Number>& couplingMatrix_mediator_BIG =
289  *m_couplingMatrixMap_mediator_BIG[micro_name];
290  libMesh::PetscMatrix<libMesh::Number>& couplingMatrix_mediator_micro =
291  *m_couplingMatrixMap_mediator_micro[micro_name];
292  libMesh::PetscMatrix<libMesh::Number>& couplingMatrix_mediator_mediator =
293  *m_couplingMatrixMap_mediator_mediator[micro_name];
294 
295  const libMesh::Parallel::Communicator& WorldComm = couplingMatrix_mediator_micro.comm();
296 
297  // Values of the coupling constants
298  double L2_coupling_const = m_coupling_rigidityMap[micro_name]
299  / (m_coupling_widthMap[micro_name]
300  * m_coupling_widthMap[micro_name]);
301  double H1_coupling_const = m_coupling_rigidityMap[micro_name];
302 
303  // DoF vectors and ranges
304  mediator_addresses.set_DoFs();
305  BIG_addresses.set_DoFs();
306  micro_addresses.set_DoFs();
307  inter_addresses.set_DoFs();
308 
309  R_BIG_addresses.set_DoFs();
310  R_micro_addresses.set_DoFs();
311 
312  // Local matrix
313  coupling_matrices_3 Me_mediator_R_micro;
314  coupling_matrices_3 Me_mediator_R_BIG;
315  coupling_matrices_3 Me_mediator_mediator;
316 
317  // If all elements are of the same type, do the index "extraction",
318  // the matrices resizes and repositions here
319  if (bSameElemsType)
320  {
321  Me_mediator_R_micro.set_matrices(mediator_addresses, R_micro_addresses);
322  Me_mediator_R_BIG.set_matrices(mediator_addresses, R_BIG_addresses);
323  Me_mediator_mediator.set_matrices(mediator_addresses,
324  mediator_addresses);
325 
326  // Restart the element
327 
328  const libMesh::Elem* elem_inter =
329  *inter_addresses.mesh.active_local_elements_begin();
330  inter_addresses.fe_unique_ptr->reinit(elem_inter);
331  n_quadrature_pts = inter_addresses.fe_unique_ptr->n_quadrature_points();
332 
333  corrected_phi_R_BIG.resize(R_BIG_addresses.n_dofs_u,
334  std::vector<libMesh::Real>(n_quadrature_pts, 0));
335  corrected_phi_R_micro.resize(R_micro_addresses.n_dofs_u,
336  std::vector<libMesh::Real>(n_quadrature_pts, 0));
337  corrected_phi_mediator.resize(mediator_addresses.n_dofs_u,
338  std::vector<libMesh::Real>(n_quadrature_pts, 0));
339 
340  if (m_bUseH1Coupling[micro_name])
341  {
342  corrected_dphi_R_BIG.resize(R_BIG_addresses.n_dofs_u,
343  std::vector<libMesh::RealGradient>(n_quadrature_pts));
344  corrected_dphi_mediator.resize(mediator_addresses.n_dofs_u,
345  std::vector<libMesh::RealGradient>(n_quadrature_pts));
346  corrected_dphi_R_micro.resize(R_micro_addresses.n_dofs_u,
347  std::vector<libMesh::RealGradient>(n_quadrature_pts));
348  }
349  }
350 
351  // Vectors containing the lambda weights
352  std::vector<std::vector<libMesh::Real> > lambda_weight_mediator(
353  mediator_addresses.n_dofs_u,
354  std::vector<libMesh::Real>(inter_addresses.n_dofs_u, 0));
355  std::vector<std::vector<libMesh::Real> > lambda_weight_R_BIG(
356  R_BIG_addresses.n_dofs_u,
357  std::vector<libMesh::Real>(inter_addresses.n_dofs_u, 0));
358  std::vector<std::vector<libMesh::Real> > lambda_weight_R_micro(
359  R_micro_addresses.n_dofs_u,
360  std::vector<libMesh::Real>(inter_addresses.n_dofs_u, 0));
361 
362  // Initialize global matrix
363  prepare_coupling_preallocation(
364  couplingMatrix_mediator_BIG,
365  mediator_addresses,
366  BIG_addresses,
367  inter_table_mediator_BIG
368  );
369 
370  prepare_coupling_preallocation(
371  couplingMatrix_mediator_micro,
372  mediator_addresses,
373  micro_addresses,
374  inter_table_mediator_micro
375  );
376 
377  couplingMatrix_mediator_mediator.attach_dof_map(mediator_addresses.dof_map);
378  couplingMatrix_mediator_mediator.init();
379 
380  // Intersection indexes and iterators
381  int inter_idx = -1;
382  int elem_restrict_BIG_idx = -1;
383  int elem_restrict_micro_idx = -1;
384 
385  int elem_BIG_idx = -1;
386  int elem_micro_idx = -1;
387  int elem_mediator_idx = -1;
388  int elem_inter_idx = -1;
389 
390  std::pair<int, int> restrict_idx_pair;
391  std::pair<int, int> idx_pair;
392 
393  // DEBUG CODE !!!
394  int DEBUG_max_inter_idx = -1;
395  int DEBUG_nb_of_inter_elems = 0;
396  double DEBUG_vol = 0;
397 
398  // For each intersection
399  libMesh::MeshBase::const_element_iterator inter_elIt =
400  inter_addresses.mesh.active_local_elements_begin();
401  const libMesh::MeshBase::const_element_iterator end_inter_elIt =
402  inter_addresses.mesh.active_local_elements_end();
403 
404  for( ; inter_elIt != end_inter_elIt; ++inter_elIt )
405  {
406  const libMesh::Elem* elem_inter = *inter_elIt;
407 
408  // Get the intersection element idx
409  elem_inter_idx = elem_inter->id();
410  inter_idx = local_intersection_meshI_to_inter_map.at(elem_inter_idx);
411 
412  restrict_idx_pair = full_intersection_restricted_pairs_map.at(inter_idx);
413  elem_restrict_BIG_idx = restrict_idx_pair.first;
414  elem_restrict_micro_idx = restrict_idx_pair.second;
415 
416  idx_pair = full_intersection_pairs_map.at(inter_idx);
417  elem_BIG_idx = idx_pair.first;
418  elem_micro_idx = idx_pair.second;
419 
420  // TODO For now, we suppose that BIG contains the mediator. So ...
421  elem_mediator_idx = elem_restrict_BIG_idx;
422 
423  const libMesh::Elem* elem_mediator =
424  mediator_addresses.mesh.elem(elem_mediator_idx);
425  const libMesh::Elem* elem_R_BIG =
426  R_BIG_addresses.mesh.elem(elem_restrict_BIG_idx);
427  const libMesh::Elem* elem_R_micro =
428  R_micro_addresses.mesh.elem(elem_restrict_micro_idx);
429 
430  mediator_addresses.set_DoFs(elem_mediator_idx);
431  R_BIG_addresses.set_DoFs(elem_restrict_BIG_idx);
432  R_micro_addresses.set_DoFs(elem_restrict_micro_idx);
433 
434  BIG_addresses.set_DoFs(elem_BIG_idx);
435  micro_addresses.set_DoFs(elem_micro_idx);
436 
437  // Resize dense matrix, if needed
438  if (!bSameElemsType)
439  {
440  Me_mediator_R_micro.set_matrices(mediator_addresses,
441  R_micro_addresses);
442  Me_mediator_R_BIG.set_matrices(mediator_addresses,
443  R_BIG_addresses);
444  Me_mediator_mediator.set_matrices(mediator_addresses,
445  mediator_addresses);
446 
447  // Set up corrected shape vectors
448  inter_addresses.fe_unique_ptr->reinit(elem_inter);
449  n_quadrature_pts =
450  inter_addresses.fe_unique_ptr->n_quadrature_points();
451 
452  corrected_phi_R_BIG.resize(R_BIG_addresses.n_dofs_u,
453  std::vector<libMesh::Real>(n_quadrature_pts, 0));
454  corrected_phi_mediator.resize(mediator_addresses.n_dofs_u,
455  std::vector<libMesh::Real>(n_quadrature_pts, 0));
456  corrected_phi_R_micro.resize(R_micro_addresses.n_dofs_u,
457  std::vector<libMesh::Real>(n_quadrature_pts, 0));
458 
459  if (m_bUseH1Coupling[micro_name])
460  {
461  corrected_dphi_R_BIG.resize(R_BIG_addresses.n_dofs_u,
462  std::vector<libMesh::RealGradient>(n_quadrature_pts));
463  corrected_dphi_mediator.resize(mediator_addresses.n_dofs_u,
464  std::vector<libMesh::RealGradient>(n_quadrature_pts));
465  corrected_dphi_R_micro.resize(R_micro_addresses.n_dofs_u,
466  std::vector<libMesh::RealGradient>(n_quadrature_pts));
467  }
468  }
469 
470  Me_mediator_R_micro.zero();
471  Me_mediator_R_BIG.zero();
472  Me_mediator_mediator.zero();
473 
474  // DEBUG!!!
475  if (MASTER_debug_coupling_assemble)
476  {
477  if (elem_inter_idx > DEBUG_max_inter_idx)
478  {
479  DEBUG_max_inter_idx = elem_inter_idx;
480  }
481  }
482 
483  // DEBUG !!!
484  if (MASTER_debug_coupling_assemble)
485  {
486  DEBUG_vol += elem_inter->volume();
487  ++DEBUG_nb_of_inter_elems;
488  }
489 
490  // Restart the element
491  inter_addresses.fe_unique_ptr->reinit(elem_inter);
492 
493  get_lambdas(R_BIG_addresses.dim, R_BIG_addresses.fe_type, elem_R_BIG,
494  quad_points_inter, quad_points_reference,
495  lambda_weight_R_BIG);
496 
497  get_lambdas(R_micro_addresses.dim, R_micro_addresses.fe_type,
498  elem_R_micro, quad_points_inter, quad_points_reference,
499  lambda_weight_R_micro);
500 
501  get_lambdas(mediator_addresses.dim, mediator_addresses.fe_type,
502  elem_mediator, quad_points_inter, quad_points_reference,
503  lambda_weight_mediator);
504 
505  set_corrected_shapes(lambda_weight_R_BIG, phi_inter,
506  corrected_phi_R_BIG);
507  set_corrected_shapes(lambda_weight_R_micro, phi_inter,
508  corrected_phi_R_micro);
509  set_corrected_shapes(lambda_weight_mediator, phi_inter,
510  corrected_phi_mediator);
511 
512  if (m_bUseH1Coupling[micro_name])
513  {
514  set_corrected_shape_gradients(lambda_weight_R_BIG, dphi_inter,
515  corrected_dphi_R_BIG);
516  set_corrected_shape_gradients(lambda_weight_R_micro, dphi_inter,
517  corrected_dphi_R_micro);
518  set_corrected_shape_gradients(lambda_weight_mediator,
519  dphi_inter, corrected_dphi_mediator);
520  }
521 
522  // For each quadrature point determinate the sub-matrices elements
523  for (unsigned int qp = 0; qp < inter_addresses.qrule.n_points();
524  qp++)
525  {
526  // Mediator -> micro coupling
527  Me_mediator_R_micro.build_L2_coupling_matrix(mediator_addresses,
528  R_micro_addresses, qp, corrected_phi_mediator,
529  corrected_phi_R_micro, JxW, L2_coupling_const);
530 
531  // Mediator -> BIG coupling
532  Me_mediator_R_BIG.build_L2_coupling_matrix(mediator_addresses,
533  R_BIG_addresses, qp, corrected_phi_mediator,
534  corrected_phi_R_BIG, JxW, L2_coupling_const);
535 
536  // Mediator -> Mediator coupling
537  Me_mediator_mediator.build_L2_coupling_matrix(
538  mediator_addresses, mediator_addresses, qp,
539  corrected_phi_mediator, corrected_phi_mediator, JxW,
540  L2_coupling_const);
541 
542  if (m_bUseH1Coupling[micro_name])
543  {
544  // Then we must also build the strain terms
545 
546  // Mediator -> micro coupling
547  Me_mediator_R_micro.add_H1_coupling_matrix(mediator_addresses,
548  R_micro_addresses, qp, corrected_dphi_mediator,
549  corrected_dphi_R_micro, JxW, H1_coupling_const);
550 
551  Me_mediator_R_BIG.add_H1_coupling_matrix(mediator_addresses,
552  R_BIG_addresses, qp, corrected_dphi_mediator,
553  corrected_dphi_R_BIG, JxW, H1_coupling_const);
554 
555  Me_mediator_mediator.add_H1_coupling_matrix(
556  mediator_addresses, mediator_addresses, qp,
557  corrected_dphi_mediator, corrected_dphi_mediator,
558  JxW, H1_coupling_const);
559  }
560  }
561 
562  couplingMatrix_mediator_micro.add_matrix(Me_mediator_R_micro.Me,
563  mediator_addresses.dof_indices, micro_addresses.dof_indices);
564  couplingMatrix_mediator_BIG.add_matrix(Me_mediator_R_BIG.Me,
565  mediator_addresses.dof_indices, BIG_addresses.dof_indices);
566  couplingMatrix_mediator_mediator.add_matrix(Me_mediator_mediator.Me,
567  mediator_addresses.dof_indices, mediator_addresses.dof_indices);
568  }
569 
570  WorldComm.barrier();
571 
572  couplingMatrix_mediator_micro.close();
573  couplingMatrix_mediator_BIG.close();
574  couplingMatrix_mediator_mediator.close();
575 
576  print_matrix_dim(couplingMatrix_mediator_micro);
577  print_matrix_dim(couplingMatrix_mediator_BIG);
578  print_matrix_dim(couplingMatrix_mediator_mediator);
579 
580  if (MASTER_debug_coupling_assemble)
581  {
582  std::cout << "> COUPLING DEBUG !!! " << std::endl;
583  std::cout << ">" << std::endl;
584  std::cout << "> DEBUG_max_inter_idx = " << DEBUG_max_inter_idx
585  << std::endl;
586  std::cout << "> DEBUG_nb_of_inter_elems = "
587  << DEBUG_nb_of_inter_elems << std::endl;
588  std::cout << "> DEBUG_vol = " << DEBUG_vol << std::endl
589  << std::endl;
590  }
591  };
592 
593  void assemble_coupling_matrices::check_coupling_construction_3D_parallel(
594  const std::string BIG_name,
595  const std::string micro_name,
596  const std::string inter_name,
597  const std::string mediator_name,
598 
599  const libMesh::MeshBase& mesh_R_BIG,
600  const libMesh::MeshBase& mesh_R_micro,
601 
602  const std::unordered_map<int,std::pair<int,int> >&
603  full_intersection_pairs_map,
604  const std::unordered_map<int,std::pair<int,int> >&
605  full_intersection_restricted_pairs_map,
606  const std::unordered_map<int,int>&
607  local_intersection_meshI_to_inter_map,
608  const std::unordered_multimap<int,int>& inter_table_mediator_BIG,
609  const std::unordered_multimap<int,int>& inter_table_mediator_micro,
610 
611  const std::string BIG_type,
612  const std::string micro_type,
613  bool bSameElemsType)
614  {
615  // TODO : make it possible to invert the algorithm's systems!
616  // Addresses to the eq. systems
617  libMesh::EquationSystems& mediator_eq_system =
618  *m_mediator_EquationSystemMap[mediator_name];
619  libMesh::EquationSystems& BIG_eq_system =
620  *m_BIG_EquationSystem.second;
621  libMesh::EquationSystems& micro_eq_system =
622  *m_micro_EquationSystemMap[micro_name];
623  libMesh::EquationSystems& inter_eq_system =
624  *m_inter_EquationSystemMap[inter_name];
625 
626  libMesh::EquationSystems& R_BIG_eq_system =
627  *m_R_BIG_EquationSystem.second;
628  libMesh::EquationSystems& R_micro_eq_system =
629  *m_R_micro_EquationSystemMap[micro_name];
630 
631  // First, test if all the systems have an elasticity model and variable set
632  homemade_assert_msg(micro_eq_system.has_system(micro_type),
633  " Micro equation systems is missing a system type!");
634  homemade_assert_msg(BIG_eq_system.has_system(BIG_type),
635  " Macro equation systems is missing a system type!");
636  homemade_assert_msg(R_micro_eq_system.has_system("Elasticity"),
637  " Restricted micro equation systems missing \"Elasticity\" system!");
638  homemade_assert_msg(R_BIG_eq_system.has_system("Elasticity"),
639  " Restricted macro equation systems missing \"Elasticity\" system!");
640  homemade_assert_msg(inter_eq_system.has_system("Elasticity"),
641  " Intersection equation systems missing \"Elasticity\" system!");
642  homemade_assert_msg(mediator_eq_system.has_system("Elasticity"),
643  " Mediatored equation systems missing \"Elasticity\" system!");
644 
645  // Systems and vars
646  libMesh::System& volume_mediator_system =
647  libMesh::cast_ref<libMesh::System&>(mediator_eq_system.get_system("Elasticity"));
648 
649  libMesh::System& volume_BIG_system =
650  libMesh::cast_ref<libMesh::System&>(BIG_eq_system.get_system<libMesh::ExplicitSystem>(BIG_type));
651 
652  libMesh::System& volume_micro_system =
653  libMesh::cast_ref<libMesh::System&>(micro_eq_system.get_system<libMesh::ExplicitSystem>(micro_type));
654 
655  libMesh::System& volume_inter_system =
656  inter_eq_system.get_system<libMesh::ExplicitSystem>("Elasticity");
657 
658  libMesh::System& volume_R_BIG_system =
659  R_BIG_eq_system.get_system<libMesh::ExplicitSystem>("Elasticity");
660 
661  libMesh::System& volume_R_micro_system =
662  R_micro_eq_system.get_system<libMesh::ExplicitSystem>("Elasticity");
663 
664  libMesh_fe_addresses_3 mediator_addresses(volume_mediator_system);
665  libMesh_fe_addresses_3 BIG_addresses(volume_BIG_system);
666  libMesh_fe_addresses_3 micro_addresses(volume_micro_system);
667  libMesh_fe_addresses_3 inter_addresses(volume_inter_system);
668 
669  libMesh_fe_addresses_3 R_BIG_addresses(volume_R_BIG_system);
670  libMesh_fe_addresses_3 R_micro_addresses(volume_R_micro_system);
671 
672  // Vector that will keep the quadrature points
673  const std::vector<libMesh::Point>& quad_points_inter =
674  inter_addresses.fe_unique_ptr->get_xyz();
675  std::vector<libMesh::Point> quad_points_reference;
676 
677  // Jacobians
678  const std::vector<libMesh::Real>& JxW =
679  inter_addresses.fe_unique_ptr->get_JxW();
680 
681  // Shape functions
682  const std::vector<std::vector<libMesh::Real> >& phi_inter =
683  inter_addresses.fe_unique_ptr->get_phi();
684  std::vector<std::vector<libMesh::Real> > corrected_phi_R_BIG;
685  std::vector<std::vector<libMesh::Real> > corrected_phi_R_micro;
686  std::vector<std::vector<libMesh::Real> > corrected_phi_mediator;
687 
688  // Shape functions gradients
689  const std::vector<std::vector<libMesh::RealGradient> >& dphi_inter =
690  inter_addresses.fe_unique_ptr->get_dphi();
691  std::vector<std::vector<libMesh::RealGradient> > corrected_dphi_R_BIG;
692  std::vector<std::vector<libMesh::RealGradient> > corrected_dphi_R_micro;
693  std::vector<std::vector<libMesh::RealGradient> > corrected_dphi_mediator;
694 
695  unsigned int n_quadrature_pts = 0;
696 
697  // Addresses to the matrices
698  libMesh::PetscMatrix<libMesh::Number> couplingMatrix_mediator_BIG(mesh_R_BIG.comm());
699  libMesh::PetscMatrix<libMesh::Number> couplingMatrix_mediator_micro(mesh_R_BIG.comm());
700  libMesh::PetscMatrix<libMesh::Number> couplingMatrix_mediator_mediator(mesh_R_BIG.comm());
701 
702  const libMesh::Parallel::Communicator& WorldComm = couplingMatrix_mediator_micro.comm();
703 
704  // DoF vectors and ranges
705  mediator_addresses.set_DoFs();
706  BIG_addresses.set_DoFs();
707  micro_addresses.set_DoFs();
708  inter_addresses.set_DoFs();
709 
710  R_BIG_addresses.set_DoFs();
711  R_micro_addresses.set_DoFs();
712 
713  // Local matrix
714  coupling_matrices_3 Me_mediator_R_micro;
715  coupling_matrices_3 Me_mediator_R_BIG;
716  coupling_matrices_3 Me_mediator_mediator;
717 
718  // If all elements are of the same type, do the index "extraction",
719  // the matrices resizes and repositions here
720  if (bSameElemsType)
721  {
722  Me_mediator_R_micro.set_matrices(mediator_addresses, R_micro_addresses);
723  Me_mediator_R_BIG.set_matrices(mediator_addresses, R_BIG_addresses);
724  Me_mediator_mediator.set_matrices(mediator_addresses,
725  mediator_addresses);
726 
727  // Restart the element
728 
729  const libMesh::Elem* elem_inter =
730  *inter_addresses.mesh.active_local_elements_begin();
731  inter_addresses.fe_unique_ptr->reinit(elem_inter);
732  n_quadrature_pts = inter_addresses.fe_unique_ptr->n_quadrature_points();
733 
734  corrected_phi_R_BIG.resize(R_BIG_addresses.n_dofs_u,
735  std::vector<libMesh::Real>(n_quadrature_pts, 0));
736  corrected_phi_R_micro.resize(R_micro_addresses.n_dofs_u,
737  std::vector<libMesh::Real>(n_quadrature_pts, 0));
738  corrected_phi_mediator.resize(mediator_addresses.n_dofs_u,
739  std::vector<libMesh::Real>(n_quadrature_pts, 0));
740 
741 
742  corrected_dphi_R_BIG.resize(R_BIG_addresses.n_dofs_u,
743  std::vector<libMesh::RealGradient>(n_quadrature_pts));
744  corrected_dphi_mediator.resize(mediator_addresses.n_dofs_u,
745  std::vector<libMesh::RealGradient>(n_quadrature_pts));
746  corrected_dphi_R_micro.resize(R_micro_addresses.n_dofs_u,
747  std::vector<libMesh::RealGradient>(n_quadrature_pts));
748 
749  }
750 
751  // Vectors containing the lambda weights
752  std::vector<std::vector<libMesh::Real> > lambda_weight_mediator(
753  mediator_addresses.n_dofs_u,
754  std::vector<libMesh::Real>(inter_addresses.n_dofs_u, 0));
755  std::vector<std::vector<libMesh::Real> > lambda_weight_R_BIG(
756  R_BIG_addresses.n_dofs_u,
757  std::vector<libMesh::Real>(inter_addresses.n_dofs_u, 0));
758  std::vector<std::vector<libMesh::Real> > lambda_weight_R_micro(
759  R_micro_addresses.n_dofs_u,
760  std::vector<libMesh::Real>(inter_addresses.n_dofs_u, 0));
761 
762  // Initialize global matrix
763  prepare_coupling_preallocation(
764  couplingMatrix_mediator_BIG,
765  mediator_addresses,
766  BIG_addresses,
767  inter_table_mediator_BIG
768  );
769 
770  prepare_coupling_preallocation(
771  couplingMatrix_mediator_micro,
772  mediator_addresses,
773  micro_addresses,
774  inter_table_mediator_micro
775  );
776 
777  couplingMatrix_mediator_mediator.attach_dof_map(mediator_addresses.dof_map);
778  couplingMatrix_mediator_mediator.init();
779 
780  std::cout << " ----------------- " << std::endl;
781  // Intersection indexes and iterators
782  int inter_idx = -1;
783  int elem_restrict_BIG_idx = -1;
784  int elem_restrict_micro_idx = -1;
785 
786  int elem_BIG_idx = -1;
787  int elem_micro_idx = -1;
788  int elem_mediator_idx = -1;
789  int elem_inter_idx = -1;
790 
791  std::pair<int, int> restrict_idx_pair;
792  std::pair<int, int> idx_pair;
793 
794  // DEBUG CODE !!!
795  int DEBUG_max_inter_idx = -1;
796  int DEBUG_nb_of_inter_elems = 0;
797  double DEBUG_vol = 0;
798 
799  // For each intersection
800  libMesh::MeshBase::const_element_iterator inter_elIt =
801  inter_addresses.mesh.active_local_elements_begin();
802  const libMesh::MeshBase::const_element_iterator end_inter_elIt =
803  inter_addresses.mesh.active_local_elements_end();
804 
805  for( ; inter_elIt != end_inter_elIt; ++inter_elIt )
806  {
807  const libMesh::Elem* elem_inter = *inter_elIt;
808 
809  // Get the intersection element idx
810  elem_inter_idx = elem_inter->id();
811  inter_idx = local_intersection_meshI_to_inter_map.at(elem_inter_idx);
812 
813  restrict_idx_pair = full_intersection_restricted_pairs_map.at(inter_idx);
814  elem_restrict_BIG_idx = restrict_idx_pair.first;
815  elem_restrict_micro_idx = restrict_idx_pair.second;
816 
817  idx_pair = full_intersection_pairs_map.at(inter_idx);
818  elem_BIG_idx = idx_pair.first;
819  elem_micro_idx = idx_pair.second;
820 
821  // TODO For now, we suppose that BIG contains the mediator. So ...
822  elem_mediator_idx = elem_restrict_BIG_idx;
823 
824  const libMesh::Elem* elem_mediator =
825  mediator_addresses.mesh.elem(elem_mediator_idx);
826  const libMesh::Elem* elem_R_BIG =
827  R_BIG_addresses.mesh.elem(elem_restrict_BIG_idx);
828  const libMesh::Elem* elem_R_micro =
829  R_micro_addresses.mesh.elem(elem_restrict_micro_idx);
830 
831  mediator_addresses.set_DoFs(elem_mediator_idx);
832  R_BIG_addresses.set_DoFs(elem_restrict_BIG_idx);
833  R_micro_addresses.set_DoFs(elem_restrict_micro_idx);
834 
835  BIG_addresses.set_DoFs(elem_BIG_idx);
836  micro_addresses.set_DoFs(elem_micro_idx);
837 
838  // Resize dense matrix, if needed
839  if (!bSameElemsType)
840  {
841  Me_mediator_R_micro.set_matrices(mediator_addresses,
842  R_micro_addresses);
843  Me_mediator_R_BIG.set_matrices(mediator_addresses,
844  R_BIG_addresses);
845  Me_mediator_mediator.set_matrices(mediator_addresses,
846  mediator_addresses);
847 
848  // Set up corrected shape vectors
849  inter_addresses.fe_unique_ptr->reinit(elem_inter);
850  n_quadrature_pts =
851  inter_addresses.fe_unique_ptr->n_quadrature_points();
852 
853  corrected_phi_R_BIG.resize(R_BIG_addresses.n_dofs_u,
854  std::vector<libMesh::Real>(n_quadrature_pts, 0));
855  corrected_phi_mediator.resize(mediator_addresses.n_dofs_u,
856  std::vector<libMesh::Real>(n_quadrature_pts, 0));
857  corrected_phi_R_micro.resize(R_micro_addresses.n_dofs_u,
858  std::vector<libMesh::Real>(n_quadrature_pts, 0));
859 
860 
861  corrected_dphi_R_BIG.resize(R_BIG_addresses.n_dofs_u,
862  std::vector<libMesh::RealGradient>(n_quadrature_pts));
863  corrected_dphi_mediator.resize(mediator_addresses.n_dofs_u,
864  std::vector<libMesh::RealGradient>(n_quadrature_pts));
865  corrected_dphi_R_micro.resize(R_micro_addresses.n_dofs_u,
866  std::vector<libMesh::RealGradient>(n_quadrature_pts));
867  }
868 
869  Me_mediator_R_micro.zero();
870  Me_mediator_R_BIG.zero();
871  Me_mediator_mediator.zero();
872 
873  // DEBUG!!!
874 
875  if (elem_inter_idx > DEBUG_max_inter_idx)
876  {
877  DEBUG_max_inter_idx = elem_inter_idx;
878  }
879 
880 
881  // DEBUG !!!
882 
883  DEBUG_vol += std::abs(elem_inter->volume());
884  ++DEBUG_nb_of_inter_elems;
885 
886 
887  // Restart the element
888  inter_addresses.fe_unique_ptr->reinit(elem_inter);
889 
890  get_lambdas(R_BIG_addresses.dim, R_BIG_addresses.fe_type, elem_R_BIG,
891  quad_points_inter, quad_points_reference,
892  lambda_weight_R_BIG);
893 
894  get_lambdas(R_micro_addresses.dim, R_micro_addresses.fe_type,
895  elem_R_micro, quad_points_inter, quad_points_reference,
896  lambda_weight_R_micro);
897 
898  get_lambdas(mediator_addresses.dim, mediator_addresses.fe_type,
899  elem_mediator, quad_points_inter, quad_points_reference,
900  lambda_weight_mediator);
901 
902  set_corrected_shapes(lambda_weight_R_BIG, phi_inter,
903  corrected_phi_R_BIG);
904  set_corrected_shapes(lambda_weight_R_micro, phi_inter,
905  corrected_phi_R_micro);
906  set_corrected_shapes(lambda_weight_mediator, phi_inter,
907  corrected_phi_mediator);
908 
909  set_corrected_shape_gradients(lambda_weight_R_BIG, dphi_inter,
910  corrected_dphi_R_BIG);
911  set_corrected_shape_gradients(lambda_weight_R_micro, dphi_inter,
912  corrected_dphi_R_micro);
913  set_corrected_shape_gradients(lambda_weight_mediator,
914  dphi_inter, corrected_dphi_mediator);
915 
916 
917  // For each quadrature point determinate the sub-matrices elements
918  for (unsigned int qp = 0; qp < inter_addresses.qrule.n_points();
919  qp++)
920  {
921  // Mediator -> micro coupling
922  Me_mediator_R_micro.build_L2_coupling_matrix(mediator_addresses,
923  R_micro_addresses, qp, corrected_phi_mediator,
924  corrected_phi_R_micro, JxW, 1);
925 
926  // Mediator -> BIG coupling
927  Me_mediator_R_BIG.build_L2_coupling_matrix(mediator_addresses,
928  R_BIG_addresses, qp, corrected_phi_mediator,
929  corrected_phi_R_BIG, JxW, 1);
930 
931  // Mediator -> Mediator coupling
932  Me_mediator_mediator.build_L2_coupling_matrix(
933  mediator_addresses, mediator_addresses, qp,
934  corrected_phi_mediator, corrected_phi_mediator, JxW,
935  1);
936 
937  // Then we must also build the strain terms
938 
939  // Mediator -> micro coupling
940  Me_mediator_R_micro.add_H1_coupling_matrix(mediator_addresses,
941  R_micro_addresses, qp, corrected_dphi_mediator,
942  corrected_dphi_R_micro, JxW, 0);
943 
944  Me_mediator_R_BIG.add_H1_coupling_matrix(mediator_addresses,
945  R_BIG_addresses, qp, corrected_dphi_mediator,
946  corrected_dphi_R_BIG, JxW, 0);
947 
948  Me_mediator_mediator.add_H1_coupling_matrix(
949  mediator_addresses, mediator_addresses, qp,
950  corrected_dphi_mediator, corrected_dphi_mediator,
951  JxW, 0);
952  }
953 
954  couplingMatrix_mediator_micro.add_matrix(Me_mediator_R_micro.Me,
955  mediator_addresses.dof_indices, micro_addresses.dof_indices);
956  couplingMatrix_mediator_BIG.add_matrix(Me_mediator_R_BIG.Me,
957  mediator_addresses.dof_indices, BIG_addresses.dof_indices);
958  couplingMatrix_mediator_mediator.add_matrix(Me_mediator_mediator.Me,
959  mediator_addresses.dof_indices, mediator_addresses.dof_indices);
960  }
961 
962  WorldComm.barrier();
963 
964  couplingMatrix_mediator_micro.close();
965  couplingMatrix_mediator_BIG.close();
966  couplingMatrix_mediator_mediator.close();
967 
968  std::cout << " ----------------- " << std::endl;
969 
970  Vec row_sum_micro, row_sum_BIG, row_sum_mediator;
971  MatCreateVecs(couplingMatrix_mediator_micro.mat(),&row_sum_micro,NULL);
972  MatCreateVecs(couplingMatrix_mediator_BIG.mat(),&row_sum_BIG,NULL);
973  MatCreateVecs(couplingMatrix_mediator_mediator.mat(),&row_sum_mediator,NULL);
974 
975  MatGetRowSum(couplingMatrix_mediator_micro.mat(),row_sum_micro);
976  MatGetRowSum(couplingMatrix_mediator_BIG.mat(),row_sum_BIG);
977  MatGetRowSum(couplingMatrix_mediator_mediator.mat(),row_sum_mediator);
978 
979  PetscScalar coupl_vol_micro = 0, coupl_vol_BIG = 0, coupl_vol_mediator = 0;
980 
981  VecSum(row_sum_micro,&coupl_vol_micro);
982  VecSum(row_sum_BIG,&coupl_vol_BIG);
983  VecSum(row_sum_mediator,&coupl_vol_mediator);
984 
985  mesh_R_BIG.comm().sum(DEBUG_vol);
986 
987  std::cout << "> COUPLING DEBUG !!! " << std::endl;
988  std::cout << ">" << std::endl;
989  std::cout << "> DEBUG_max_inter_idx = " << DEBUG_max_inter_idx
990  << std::endl;
991  std::cout << "> DEBUG_nb_of_inter_elems = "
992  << DEBUG_nb_of_inter_elems << std::endl;
993  std::cout << "> DEBUG_mesh_vol = " << DEBUG_vol
994  << std::endl;
995  std::cout << "> DEBUG_micro_vol = " << coupl_vol_micro << std::endl;
996  std::cout << "> DEBUG_big_vol = " << coupl_vol_BIG << std::endl;
997  std::cout << "> DEBUG_mediator_vol = " << coupl_vol_mediator << std::endl
998  << std::endl;
999 
1000  VecDestroy(&row_sum_micro);
1001  VecDestroy(&row_sum_BIG);
1002  VecDestroy(&row_sum_mediator);
1003  };
1004 };
carl::libMesh_fe_addresses_3::fe_type
libMesh::FEType fe_type
Definition: assemble_coupling.h:62
carl::assemble_coupling_matrices::check_coupling_construction_3D_parallel
void check_coupling_construction_3D_parallel(const std::string BIG_name, const std::string micro_name, const std::string inter_name, const std::string mediator_name, const libMesh::MeshBase &mesh_R_BIG, const libMesh::MeshBase &mesh_R_micro, const std::unordered_map< int, std::pair< int, int > > &full_intersection_pairs_map, const std::unordered_map< int, std::pair< int, int > > &full_intersection_restricted_pairs_map, const std::unordered_map< int, int > &local_intersection_meshI_to_inter_map, const std::unordered_multimap< int, int > &inter_table_mediator_BIG, const std::unordered_multimap< int, int > &inter_table_mediator_micro, const std::string BIG_type="Elasticity", const std::string micro_type="Elasticity", bool bSameElemsType=true)
Definition: assemble_functions_coupling.cpp:593
carl::coupling_matrices_3::add_H1_coupling_matrix
void add_H1_coupling_matrix(const libMesh_fe_addresses_3 &system_type_AAA, const libMesh_fe_addresses_3 &system_type_BBB, int qp, const std::vector< std::vector< libMesh::RealGradient > > &corrected_dphi_sysAAA, const std::vector< std::vector< libMesh::RealGradient > > &corrected_dphi_sysBBB, const std::vector< libMesh::Real > &JxW, const libMesh::Number H1_coupling_const)
Definition: assemble_coupling.h:176
carl::libMesh_fe_addresses_3::n_dofs
unsigned int n_dofs
Definition: assemble_coupling.h:73
carl::assemble_coupling_matrices::get_lambdas
void get_lambdas(const unsigned int dim, const libMesh::FEType &fe_t, const libMesh::Elem *base_elem, const std::vector< libMesh::Point > &phys_points, std::vector< libMesh::Point > &ref_points, std::vector< std::vector< libMesh::Real > > &lambda_weights)
Definition: assemble_coupling.cpp:291
carl::assemble_coupling_matrices::m_mediator_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_mediator_EquationSystemMap
Definition: assemble_coupling.h:262
MASTER_debug_coupling_assemble
const bool MASTER_debug_coupling_assemble
Definition: assemble_coupling.h:18
carl::assemble_coupling_matrices::m_BIG_EquationSystem
std::pair< std::string, libMesh::EquationSystems * > m_BIG_EquationSystem
Definition: assemble_coupling.h:255
carl::coupling_matrices_3::set_matrices
void set_matrices(libMesh_fe_addresses_3 &system_type_AAA, libMesh_fe_addresses_3 &system_type_BBB)
Definition: assemble_coupling.h:115
carl::assemble_coupling_matrices::m_coupling_widthMap
std::map< std::string, double > m_coupling_widthMap
Definition: assemble_coupling.h:278
carl::libMesh_fe_addresses_3::dof_map
const libMesh::DofMap & dof_map
Definition: assemble_coupling.h:61
carl::assemble_coupling_matrices::m_coupling_rigidityMap
std::map< std::string, double > m_coupling_rigidityMap
Definition: assemble_coupling.h:277
carl::assemble_coupling_matrices::m_inter_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_inter_EquationSystemMap
Definition: assemble_coupling.h:261
carl::assemble_coupling_matrices::m_R_micro_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_R_micro_EquationSystemMap
Definition: assemble_coupling.h:259
carl::assemble_coupling_matrices::set_corrected_shape_gradients
void set_corrected_shape_gradients(const std::vector< std::vector< libMesh::Real > > &lambda_weights, const std::vector< std::vector< libMesh::RealGradient > > &phi_inter, std::vector< std::vector< libMesh::RealGradient > > &phi_corrected)
Definition: assemble_coupling.cpp:250
assemble_coupling.h
carl::coupling_matrices_3
Definition: assemble_coupling.h:93
carl::assemble_coupling_matrices::m_micro_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_micro_EquationSystemMap
Definition: assemble_coupling.h:258
carl::libMesh_fe_addresses_3::n_dofs_u
unsigned int n_dofs_u
Definition: assemble_coupling.h:74
carl::coupling_matrices_3::zero
void zero()
Definition: assemble_coupling.h:240
carl::assemble_coupling_matrices::assemble_coupling_elasticity_3D_parallel
void assemble_coupling_elasticity_3D_parallel(const std::string BIG_name, const std::string micro_name, const std::string inter_name, const std::string mediator_name, const libMesh::MeshBase &mesh_R_BIG, const libMesh::MeshBase &mesh_R_micro, const std::unordered_map< int, std::pair< int, int > > &full_intersection_pairs_map, const std::unordered_map< int, std::pair< int, int > > &full_intersection_restricted_pairs_map, const std::unordered_map< int, int > &local_intersection_meshI_to_inter_map, const std::unordered_multimap< int, int > &inter_table_mediator_BIG, const std::unordered_multimap< int, int > &inter_table_mediator_micro, const std::string BIG_type="Elasticity", const std::string micro_type="Elasticity", bool bSameElemsType=true)
Definition: assemble_functions_coupling.cpp:183
carl::libMesh_fe_addresses_3::mesh
const libMesh::MeshBase & mesh
Definition: assemble_coupling.h:56
carl::libMesh_fe_addresses_3::dim
const unsigned int dim
Definition: assemble_coupling.h:57
carl::assemble_coupling_matrices::m_R_BIG_EquationSystem
std::pair< std::string, libMesh::EquationSystems * > m_R_BIG_EquationSystem
Definition: assemble_coupling.h:256
carl::coupling_matrices_3::Me
libMesh::DenseMatrix< libMesh::Number > Me
Definition: assemble_coupling.h:96
carl::assemble_coupling_matrices::m_bUseH1Coupling
std::map< std::string, bool > m_bUseH1Coupling
Definition: assemble_coupling.h:270
carl::coupling_matrices_3::build_L2_coupling_matrix
void build_L2_coupling_matrix(const libMesh_fe_addresses_3 &system_type_AAA, const libMesh_fe_addresses_3 &system_type_BBB, int qp, const std::vector< std::vector< libMesh::Real > > &corrected_phi_AAA, const std::vector< std::vector< libMesh::Real > > &corrected_phi_BBB, const std::vector< libMesh::Real > &JxW, double L2_coupling_const)
Definition: assemble_coupling.h:157
carl::libMesh_fe_addresses_3::dof_indices
std::vector< libMesh::dof_id_type > dof_indices
Definition: assemble_coupling.h:68
carl::libMesh_fe_addresses_3::set_DoFs
void set_DoFs(int idx=0)
Definition: assemble_coupling.h:78
carl::assemble_coupling_matrices::m_couplingMatrixMap_mediator_mediator
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_mediator
Definition: assemble_coupling.h:267
carl::assemble_coupling_matrices::m_couplingMatrixMap_mediator_BIG
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_BIG
Definition: assemble_coupling.h:266
carl::assemble_coupling_matrices::set_corrected_shapes
void set_corrected_shapes(const std::vector< std::vector< libMesh::Real > > &lambda_weights, const std::vector< std::vector< libMesh::Real > > &phi_inter, std::vector< std::vector< libMesh::Real > > &phi_corrected)
Definition: assemble_coupling.cpp:209
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
carl::assemble_coupling_matrices::m_couplingMatrixMap_mediator_micro
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_micro
Definition: assemble_coupling.h:265
carl::print_matrix_dim
void print_matrix_dim(libMesh::PetscMatrix< libMesh::Number > &CouplingTestMatrix, bool bDetailed=false)
Definition: PETSC_matrix_operations.cpp:165
carl::assemble_coupling_matrices::prepare_coupling_preallocation
void prepare_coupling_preallocation(libMesh::PetscMatrix< libMesh::Number > &coupling_matrix, libMesh_fe_addresses_3 &row_addresses, libMesh_fe_addresses_3 &col_addresses, const std::unordered_multimap< int, int > &inter_table)
Definition: assemble_functions_coupling.cpp:6
carl::libMesh_fe_addresses_3::qrule
libMesh::QGauss qrule
Definition: assemble_coupling.h:64
carl::MPI_reduce_vector
void MPI_reduce_vector(std::vector< T > &r, int root, const libMesh::Parallel::Communicator &Comm)
Definition: mpi_carl_tools.h:22
carl::libMesh_fe_addresses_3
Definition: assemble_coupling.h:23
carl::libMesh_fe_addresses_3::fe_unique_ptr
libMesh::UniquePtr< libMesh::FEBase > fe_unique_ptr
Definition: assemble_coupling.h:63