CArl
Code Arlequin / C++ implementation
Public Member Functions | Protected Types | Protected Attributes | List of all members
carl::assemble_coupling_matrices Class Reference

#include <assemble_coupling.h>

Public Member Functions

 assemble_coupling_matrices (const libMesh::Parallel::Communicator &comm)
 
 ~assemble_coupling_matrices ()
 
void clear ()
 
libMesh::EquationSystems & set_BIG_EquationSystem (const std::string &name, libMesh::MeshBase &BIGMesh)
 
libMesh::EquationSystems & set_Restricted_BIG_EquationSystem (const std::string &name, libMesh::MeshBase &R_BIGMesh)
 
template<typename libMesh_MatrixType >
libMesh::EquationSystems & add_micro_EquationSystem (const std::string &name, libMesh::MeshBase &microMesh)
 
libMesh::EquationSystems & add_Restricted_micro_EquationSystem (const std::string &name, libMesh::MeshBase &R_microMesh)
 
libMesh::EquationSystems & add_inter_EquationSystem (const std::string &name, libMesh::MeshBase &interMesh)
 
libMesh::EquationSystems & add_mediator_EquationSystem (const std::string &name, libMesh::MeshBase &mediatorMesh)
 
void set_coupling_parameters (const std::string &name, double coupling_rigidity, double coupling_width)
 Set physical parameters. More...
 
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)
 
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)
 
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)
 
libMesh::PetscMatrix< libMesh::Number > & get_micro_coupling_matrix (const std::string &name)
 
libMesh::PetscMatrix< libMesh::Number > & get_BIG_coupling_matrix (const std::string &name)
 
libMesh::PetscMatrix< libMesh::Number > & get_mediator_coupling_matrix (const std::string &name)
 
void print_matrix_micro_info (const std::string &name)
 
void print_matrix_BIG_info (const std::string &name)
 
void print_matrix_mediator_info (const std::string &name)
 
void print_matrices_matlab (const std::string &name, const std::string &outputRoot="coupling")
 
void print_PETSC_matrices (const std::string &name, const std::string &outputRoot="coupling")
 
void use_H1_coupling (std::string name)
 
void use_L2_coupling (std::string name)
 
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)
 
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)
 
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)
 

Protected Types

typedef std::map< std::string, libMesh::EquationSystems * >::iterator EqSystem_iterator
 
typedef std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * >::iterator Matrix_iterator
 
typedef std::map< std::string, libMesh::PetscVector< libMesh::Number > * >::iterator Vector_iterator
 

Protected Attributes

const libMesh::Parallel::Communicator * m_comm
 
std::pair< std::string, libMesh::EquationSystems * > m_BIG_EquationSystem
 
std::pair< std::string, libMesh::EquationSystems * > m_R_BIG_EquationSystem
 
std::map< std::string, libMesh::EquationSystems * > m_micro_EquationSystemMap
 
std::map< std::string, libMesh::EquationSystems * > m_R_micro_EquationSystemMap
 
std::map< std::string, libMesh::EquationSystems * > m_inter_EquationSystemMap
 
std::map< std::string, libMesh::EquationSystems * > m_mediator_EquationSystemMap
 
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_micro
 
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_BIG
 
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_mediator
 
std::map< std::string, bool > m_bUseH1Coupling
 
std::map< std::string, bool > m_bHasAssembled_micro
 
bool m_bHasAssembled_BIG
 
bool m_bUseNullSpace_BIG
 
bool m_bHasDefinedMeshRestrictions
 
std::map< std::string, double > m_coupling_rigidityMap
 
std::map< std::string, double > m_coupling_widthMap
 

Detailed Description

Definition at line 246 of file assemble_coupling.h.

Member Typedef Documentation

typedef std::map<std::string, libMesh::EquationSystems*>::iterator carl::assemble_coupling_matrices::EqSystem_iterator
protected

Definition at line 281 of file assemble_coupling.h.

typedef std::map<std::string, libMesh::PetscMatrix<libMesh::Number>*>::iterator carl::assemble_coupling_matrices::Matrix_iterator
protected

Definition at line 282 of file assemble_coupling.h.

typedef std::map<std::string, libMesh::PetscVector<libMesh::Number>*>::iterator carl::assemble_coupling_matrices::Vector_iterator
protected

Definition at line 283 of file assemble_coupling.h.

Constructor & Destructor Documentation

carl::assemble_coupling_matrices::assemble_coupling_matrices ( const libMesh::Parallel::Communicator &  comm)
inline

Definition at line 292 of file assemble_coupling.h.

292  :
293  m_comm { &comm },
294  m_bHasAssembled_BIG { false },
295  m_bHasDefinedMeshRestrictions { false }
296  {
297  };
carl::assemble_coupling_matrices::m_comm
const libMesh::Parallel::Communicator * m_comm
Definition: assemble_coupling.h:252
carl::assemble_coupling_matrices::m_bHasAssembled_BIG
bool m_bHasAssembled_BIG
Definition: assemble_coupling.h:272
carl::assemble_coupling_matrices::m_bHasDefinedMeshRestrictions
bool m_bHasDefinedMeshRestrictions
Definition: assemble_coupling.h:274
carl::assemble_coupling_matrices::~assemble_coupling_matrices ( )
inline

Definition at line 300 of file assemble_coupling.h.

301  {
302  this->clear();
303  }
carl::assemble_coupling_matrices::clear
void clear()
Definition: assemble_coupling.cpp:6

Member Function Documentation

libMesh::EquationSystems & carl::assemble_coupling_matrices::add_inter_EquationSystem ( const std::string &  name,
libMesh::MeshBase &  interMesh 
)

Definition at line 157 of file assemble_coupling.cpp.

159  {
160  libMesh::EquationSystems* EqSystemPtr = NULL;
161 
162  if (!m_inter_EquationSystemMap.count(name))
163  {
164  // Then add a new system
165  EqSystemPtr = new libMesh::EquationSystems(interMesh);
166  m_inter_EquationSystemMap.insert(std::make_pair(name, EqSystemPtr));
167  }
168  else
169  {
170  // System already exits, return the system pointer
171  std::cerr << " *** Warning: inter system " << name
172  << " already exists!" << std::endl;
173  EqSystemPtr = m_inter_EquationSystemMap[name];
174  }
175 
176  return *EqSystemPtr;
177  };
carl::assemble_coupling_matrices::m_inter_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_inter_EquationSystemMap
Definition: assemble_coupling.h:261
libMesh::EquationSystems & carl::assemble_coupling_matrices::add_mediator_EquationSystem ( const std::string &  name,
libMesh::MeshBase &  mediatorMesh 
)

Definition at line 179 of file assemble_coupling.cpp.

181  {
182  libMesh::EquationSystems* EqSystemPtr = NULL;
183 
184  if (!m_mediator_EquationSystemMap.count(name))
185  {
186  // Then add a new system
187  EqSystemPtr = new libMesh::EquationSystems(mediatorMesh);
188  m_mediator_EquationSystemMap.insert(
189  std::make_pair(name, EqSystemPtr));
190  }
191  else
192  {
193  // System already exits, return the system pointer
194  std::cerr << " *** Warning: mediator system " << name
195  << " already exists!" << std::endl;
196  EqSystemPtr = m_mediator_EquationSystemMap[name];
197  }
198 
199  return *EqSystemPtr;
200  };
carl::assemble_coupling_matrices::m_mediator_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_mediator_EquationSystemMap
Definition: assemble_coupling.h:262
template<typename libMesh_MatrixType >
libMesh::EquationSystems& carl::assemble_coupling_matrices::add_micro_EquationSystem ( const std::string &  name,
libMesh::MeshBase &  microMesh 
)
inline

Definition at line 316 of file assemble_coupling.h.

318  {
319  libMesh::EquationSystems* EqSystemPtr = NULL;
320  libMesh::PetscMatrix<libMesh::Number>* Matrix_mediator_micro_Ptr = NULL;
321  libMesh::PetscMatrix<libMesh::Number>* Matrix_mediator_BIG_Ptr = NULL;
322  libMesh::PetscMatrix<libMesh::Number>* Matrix_mediator_mediator_Ptr =
323  NULL;
324 
325  if (!m_micro_EquationSystemMap.count(name))
326  {
327  // Then add a new system
328  EqSystemPtr = new libMesh::EquationSystems(microMesh);
329 
330  Matrix_mediator_micro_Ptr = new libMesh_MatrixType(
331  microMesh.comm());
332  Matrix_mediator_BIG_Ptr = new libMesh_MatrixType(microMesh.comm());
333  Matrix_mediator_mediator_Ptr = new libMesh_MatrixType(
334  microMesh.comm());
335 
336  m_micro_EquationSystemMap.insert(std::make_pair(name, EqSystemPtr));
337  m_couplingMatrixMap_mediator_micro.insert(
338  std::make_pair(name, Matrix_mediator_micro_Ptr));
339 
340  m_couplingMatrixMap_mediator_BIG.insert(
341  std::make_pair(name, Matrix_mediator_BIG_Ptr));
342 
343  m_couplingMatrixMap_mediator_mediator.insert(
344  std::make_pair(name, Matrix_mediator_mediator_Ptr));
345 
346  m_bHasAssembled_micro.insert(std::make_pair(name, false));
347  m_bUseH1Coupling.insert(std::make_pair(name, false));
348  }
349  else
350  {
351  // System already exits, return the system pointer
352  std::cerr << " *** Warning: micro system " << name
353  << " already exists!" << std::endl;
354  EqSystemPtr = m_micro_EquationSystemMap[name];
355  }
356 
357  return *EqSystemPtr;
358  }
carl::assemble_coupling_matrices::m_micro_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_micro_EquationSystemMap
Definition: assemble_coupling.h:258
carl::assemble_coupling_matrices::m_bHasAssembled_micro
std::map< std::string, bool > m_bHasAssembled_micro
Definition: assemble_coupling.h:271
carl::assemble_coupling_matrices::m_bUseH1Coupling
std::map< std::string, bool > m_bUseH1Coupling
Definition: assemble_coupling.h:270
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::m_couplingMatrixMap_mediator_micro
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_micro
Definition: assemble_coupling.h:265
libMesh::EquationSystems & carl::assemble_coupling_matrices::add_Restricted_micro_EquationSystem ( const std::string &  name,
libMesh::MeshBase &  R_microMesh 
)

Definition at line 132 of file assemble_coupling.cpp.

134  {
135  m_bHasDefinedMeshRestrictions = true;
136 
137  libMesh::EquationSystems* EqSystemPtr = NULL;
138 
139  if (!m_R_micro_EquationSystemMap.count(name))
140  {
141  // Then add a new system
142  EqSystemPtr = new libMesh::EquationSystems(R_microMesh);
143 
144  m_R_micro_EquationSystemMap.insert(std::make_pair(name, EqSystemPtr));
145  }
146  else
147  {
148  // System already exits, return the system pointer
149  std::cerr << " *** Warning: restricted micro system " << name
150  << " already exists!" << std::endl;
151  EqSystemPtr = m_micro_EquationSystemMap[name];
152  }
153 
154  return *EqSystemPtr;
155  };
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::m_micro_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_micro_EquationSystemMap
Definition: assemble_coupling.h:258
carl::assemble_coupling_matrices::m_bHasDefinedMeshRestrictions
bool m_bHasDefinedMeshRestrictions
Definition: assemble_coupling.h:274
void carl::assemble_coupling_matrices::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 at line 183 of file assemble_functions_coupling.cpp.

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  };
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::assemble_coupling_matrices::m_coupling_widthMap
std::map< std::string, double > m_coupling_widthMap
Definition: assemble_coupling.h:278
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
carl::assemble_coupling_matrices::m_micro_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_micro_EquationSystemMap
Definition: assemble_coupling.h:258
carl::assemble_coupling_matrices::m_R_BIG_EquationSystem
std::pair< std::string, libMesh::EquationSystems * > m_R_BIG_EquationSystem
Definition: assemble_coupling.h:256
carl::assemble_coupling_matrices::m_bUseH1Coupling
std::map< std::string, bool > m_bUseH1Coupling
Definition: assemble_coupling.h:270
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
void carl::assemble_coupling_matrices::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 at line 593 of file assemble_functions_coupling.cpp.

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  };
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
carl::assemble_coupling_matrices::m_BIG_EquationSystem
std::pair< std::string, libMesh::EquationSystems * > m_BIG_EquationSystem
Definition: assemble_coupling.h:255
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
carl::assemble_coupling_matrices::m_micro_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_micro_EquationSystemMap
Definition: assemble_coupling.h:258
carl::assemble_coupling_matrices::m_R_BIG_EquationSystem
std::pair< std::string, libMesh::EquationSystems * > m_R_BIG_EquationSystem
Definition: assemble_coupling.h:256
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::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
void carl::assemble_coupling_matrices::clear ( )

Definition at line 6 of file assemble_coupling.cpp.

7  {
8  // Clean up all systems
9  libMesh::EquationSystems *EqBIGSys = m_BIG_EquationSystem.second;
10  EqBIGSys->clear();
11  delete EqBIGSys;
12  EqBIGSys = NULL;
13  m_BIG_EquationSystem.second = NULL;
14 
15  while(!m_micro_EquationSystemMap.empty())
16  {
17  EqSystem_iterator toClean = m_micro_EquationSystemMap.begin();
18 
19  libMesh::EquationSystems *EqSys = toClean->second;
20  EqSys->clear();
21  delete EqSys;
22  EqSys = NULL;
23 
24  m_micro_EquationSystemMap.erase(toClean);
25  }
26 
27  if(m_bHasDefinedMeshRestrictions)
28  {
29  libMesh::EquationSystems *EqRBIGSys = m_R_BIG_EquationSystem.second;
30  EqRBIGSys->clear();
31  delete EqRBIGSys;
32  EqRBIGSys = NULL;
33  m_R_BIG_EquationSystem.second = NULL;
34 
35  while(!m_R_micro_EquationSystemMap.empty())
36  {
37  EqSystem_iterator toClean = m_R_micro_EquationSystemMap.begin();
38 
39  libMesh::EquationSystems *EqSys = toClean->second;
40  EqSys->clear();
41  delete EqSys;
42  EqSys = NULL;
43 
44  m_R_micro_EquationSystemMap.erase(toClean);
45  }
46  }
47 
48  while(!m_inter_EquationSystemMap.empty())
49  {
50  EqSystem_iterator toClean = m_inter_EquationSystemMap.begin();
51 
52  libMesh::EquationSystems *EqSys = toClean->second;
53  EqSys->clear();
54  delete EqSys;
55  EqSys = NULL;
56 
57  m_inter_EquationSystemMap.erase(toClean);
58  }
59 
60  while(!m_mediator_EquationSystemMap.empty())
61  {
62  EqSystem_iterator toClean = m_mediator_EquationSystemMap.begin();
63 
64  libMesh::EquationSystems *EqSys = toClean->second;
65  EqSys->clear();
66  delete EqSys;
67  EqSys = NULL;
68 
69  m_mediator_EquationSystemMap.erase(toClean);
70  }
71 
72  while(!m_couplingMatrixMap_mediator_micro.empty())
73  {
74  Matrix_iterator toClean = m_couplingMatrixMap_mediator_micro.begin();
75 
76  libMesh::PetscMatrix<libMesh::Number> *Mat = toClean->second;
77  Mat->clear();
78  delete Mat;
79  Mat = NULL;
80 
81  m_couplingMatrixMap_mediator_micro.erase(toClean);
82  }
83 
84  while(!m_couplingMatrixMap_mediator_BIG.empty())
85  {
86  Matrix_iterator toClean = m_couplingMatrixMap_mediator_BIG.begin();
87 
88  libMesh::PetscMatrix<libMesh::Number> *Mat = toClean->second;
89  Mat->clear();
90  delete Mat;
91  Mat = NULL;
92 
93  m_couplingMatrixMap_mediator_BIG.erase(toClean);
94  }
95 
96  while(!m_couplingMatrixMap_mediator_mediator.empty())
97  {
98  Matrix_iterator toClean = m_couplingMatrixMap_mediator_mediator.begin();
99 
100  libMesh::PetscMatrix<libMesh::Number> *Mat = toClean->second;
101  Mat->clear();
102  delete Mat;
103  Mat = NULL;
104 
105  m_couplingMatrixMap_mediator_mediator.erase(toClean);
106  }
107  };
carl::assemble_coupling_matrices::m_mediator_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_mediator_EquationSystemMap
Definition: assemble_coupling.h:262
carl::assemble_coupling_matrices::m_BIG_EquationSystem
std::pair< std::string, libMesh::EquationSystems * > m_BIG_EquationSystem
Definition: assemble_coupling.h:255
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::m_micro_EquationSystemMap
std::map< std::string, libMesh::EquationSystems * > m_micro_EquationSystemMap
Definition: assemble_coupling.h:258
carl::assemble_coupling_matrices::m_bHasDefinedMeshRestrictions
bool m_bHasDefinedMeshRestrictions
Definition: assemble_coupling.h:274
carl::assemble_coupling_matrices::EqSystem_iterator
std::map< std::string, libMesh::EquationSystems * >::iterator EqSystem_iterator
Definition: assemble_coupling.h:281
carl::assemble_coupling_matrices::Matrix_iterator
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * >::iterator Matrix_iterator
Definition: assemble_coupling.h:282
carl::assemble_coupling_matrices::m_R_BIG_EquationSystem
std::pair< std::string, libMesh::EquationSystems * > m_R_BIG_EquationSystem
Definition: assemble_coupling.h:256
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::m_couplingMatrixMap_mediator_micro
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_micro
Definition: assemble_coupling.h:265
libMesh::PetscMatrix< libMesh::Number > & carl::assemble_coupling_matrices::get_BIG_coupling_matrix ( const std::string &  name)

Definition at line 341 of file assemble_coupling.cpp.

342  {
343  return * m_couplingMatrixMap_mediator_BIG[name];
344  }
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
void carl::assemble_coupling_matrices::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 at line 291 of file assemble_coupling.cpp.

297  {
298  // Test if we are using the correct type!
299  homemade_assert_msg(base_elem->type() == libMesh::TET4 || base_elem->type() == libMesh::HEX8,
300  " Only implemented for TET4 or HEX8 elements!");
301 
302  homemade_assert_msg(dim == 3,
303  " Only implemented for 3D!");
304 
305  unsigned int lambda_inter_size = lambda_weights[0].size();
306  unsigned int lambda_base_size = lambda_weights.size();
307 
308  // Set vectors to have same sizes
309  if(phys_points.size() != ref_points.size())
310  {
311  ref_points.resize( phys_points.size() );
312  }
313 
314  // Convert the points
315  libMesh::FEInterface::inverse_map(dim, fe_t, base_elem, phys_points, ref_points);
316 
317 
318 
319  // Calculate the lambdas
320  // -> lines : DoF from base
321  // -> cols : DoF from inter
322  for(unsigned int iii = 0; iii < lambda_base_size; ++iii)
323  {
324  for(unsigned int jjj = 0; jjj < lambda_inter_size; ++jjj)
325  {
326  lambda_weights[iii][jjj] =
327  libMesh::FE<3,libMesh::LAGRANGE>::shape(
328  base_elem->type(),
329  libMesh::FIRST,
330  iii,
331  ref_points[jjj]);
332  }
333  }
334  };
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
libMesh::PetscMatrix< libMesh::Number > & carl::assemble_coupling_matrices::get_mediator_coupling_matrix ( const std::string &  name)

Definition at line 346 of file assemble_coupling.cpp.

347  {
348  return * m_couplingMatrixMap_mediator_mediator[name];
349  }
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
libMesh::PetscMatrix< libMesh::Number > & carl::assemble_coupling_matrices::get_micro_coupling_matrix ( const std::string &  name)

Definition at line 336 of file assemble_coupling.cpp.

337  {
338  return * m_couplingMatrixMap_mediator_micro[name];
339  }
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
void carl::assemble_coupling_matrices::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 at line 6 of file assemble_functions_coupling.cpp.

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 }
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
void carl::assemble_coupling_matrices::print_matrices_matlab ( const std::string &  name,
const std::string &  outputRoot = "coupling" 
)

Definition at line 375 of file assemble_coupling.cpp.

376  {
377  libMesh::PetscMatrix<libMesh::Number>& CouplingTestMatrix_BIG =
378  * m_couplingMatrixMap_mediator_BIG[name];
379  libMesh::PetscMatrix<libMesh::Number>& CouplingTestMatrix_micro=
380  * m_couplingMatrixMap_mediator_micro[name];
381  libMesh::PetscMatrix<libMesh::Number>& CouplingTestMatrix_mediator=
382  * m_couplingMatrixMap_mediator_mediator[name];
383 
384 // Print MatLab debugging output? Variable defined at "carl_headers.h"
385 #ifdef PRINT_MATLAB_DEBUG
386  CouplingTestMatrix_BIG.print_matlab(outputRoot + "/coupling_matrix_macro.m");
387  CouplingTestMatrix_micro.print_matlab(outputRoot + "/coupling_matrix_micro.m");
388  CouplingTestMatrix_mediator.print_matlab(outputRoot + "/coupling_matrix_mediator.m");
389 #endif
390  }
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::m_couplingMatrixMap_mediator_micro
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_micro
Definition: assemble_coupling.h:265
void carl::assemble_coupling_matrices::print_matrix_BIG_info ( const std::string &  name)

Definition at line 359 of file assemble_coupling.cpp.

360  {
361  libMesh::PetscMatrix<libMesh::Number>& CouplingTestMatrix =
362  * m_couplingMatrixMap_mediator_BIG[name];
363  std::cout << "| Restrict - Macro matrix -> " << name << std::endl;
364  print_matrix(CouplingTestMatrix);
365  }
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::print_matrix
void print_matrix(libMesh::PetscMatrix< libMesh::Number > &CouplingTestMatrix)
Definition: PETSC_matrix_operations.cpp:104
void carl::assemble_coupling_matrices::print_matrix_mediator_info ( const std::string &  name)

Definition at line 367 of file assemble_coupling.cpp.

368  {
369  libMesh::PetscMatrix<libMesh::Number>& CouplingTestMatrix =
370  * m_couplingMatrixMap_mediator_mediator[name];
371  std::cout << "| Restrict - Restrict matrix -> " << name << std::endl;
372  print_matrix(CouplingTestMatrix);
373  }
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::print_matrix
void print_matrix(libMesh::PetscMatrix< libMesh::Number > &CouplingTestMatrix)
Definition: PETSC_matrix_operations.cpp:104
void carl::assemble_coupling_matrices::print_matrix_micro_info ( const std::string &  name)

Definition at line 351 of file assemble_coupling.cpp.

352  {
353  libMesh::PetscMatrix<libMesh::Number>& CouplingTestMatrix =
354  * m_couplingMatrixMap_mediator_micro[name];
355  std::cout << "| Restrict - Micro matrix -> " << name << std::endl;
356  print_matrix(CouplingTestMatrix);
357  }
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
void print_matrix(libMesh::PetscMatrix< libMesh::Number > &CouplingTestMatrix)
Definition: PETSC_matrix_operations.cpp:104
void carl::assemble_coupling_matrices::print_PETSC_matrices ( const std::string &  name,
const std::string &  outputRoot = "coupling" 
)

Definition at line 392 of file assemble_coupling.cpp.

393  {
394  libMesh::PetscMatrix<libMesh::Number>& CouplingTestMatrix_BIG =
395  * m_couplingMatrixMap_mediator_BIG[name];
396  libMesh::PetscMatrix<libMesh::Number>& CouplingTestMatrix_micro=
397  * m_couplingMatrixMap_mediator_micro[name];
398  libMesh::PetscMatrix<libMesh::Number>& CouplingTestMatrix_mediator=
399  * m_couplingMatrixMap_mediator_mediator[name];
400 
401  write_PETSC_matrix(CouplingTestMatrix_BIG.mat(), outputRoot + "/coupling_matrix_macro.petscmat",m_comm->rank(),m_comm->get());
402  write_PETSC_matrix(CouplingTestMatrix_micro.mat(), outputRoot + "/coupling_matrix_micro.petscmat",m_comm->rank(),m_comm->get());
403  write_PETSC_matrix(CouplingTestMatrix_mediator.mat(), outputRoot + "/coupling_matrix_mediator.petscmat",m_comm->rank(),m_comm->get());
404  }
carl::assemble_coupling_matrices::m_comm
const libMesh::Parallel::Communicator * m_comm
Definition: assemble_coupling.h:252
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::m_couplingMatrixMap_mediator_micro
std::map< std::string, libMesh::PetscMatrix< libMesh::Number > * > m_couplingMatrixMap_mediator_micro
Definition: assemble_coupling.h:265
carl::write_PETSC_matrix
void write_PETSC_matrix(Mat input_mat, const std::string &filename, int rank, MPI_Comm comm=PETSC_COMM_WORLD, int dim=1)
Definition: PETSC_matrix_operations.cpp:415
libMesh::EquationSystems & carl::assemble_coupling_matrices::set_BIG_EquationSystem ( const std::string &  name,
libMesh::MeshBase &  BIGMesh 
)

Definition at line 109 of file assemble_coupling.cpp.

110  {
111  libMesh::EquationSystems* EqSystemPtr = new libMesh::EquationSystems(
112  BIGMesh);
113 
114  m_BIG_EquationSystem.first = name;
115  m_BIG_EquationSystem.second = EqSystemPtr;
116 
117  return *EqSystemPtr;
118  };
carl::assemble_coupling_matrices::m_BIG_EquationSystem
std::pair< std::string, libMesh::EquationSystems * > m_BIG_EquationSystem
Definition: assemble_coupling.h:255
void carl::assemble_coupling_matrices::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 at line 250 of file assemble_coupling.cpp.

253  {
254  unsigned int lambda_corr_size = lambda_weights.size();
255  unsigned int lambda_inter_size = lambda_weights[0].size();
256 
257  unsigned int phi_inter_dof = dphi_inter.size();
258  unsigned int phi_inter_qp = dphi_inter[0].size();
259 
260  unsigned int phi_corrected_dof = dphi_corrected.size();
261  unsigned int phi_corrected_qp = dphi_corrected[0].size();
262 
263  /*
264  * lambda : [ n_dofs_corr ] x [ n_dofs_inter ]
265  * inter : [ n_dofs_inter ] x [ n_qp ]
266  * corr : [ n_dofs_corr ] x [ n_qp ]
267  */
268 
269  homemade_assert_msg(phi_inter_qp == phi_corrected_qp,
270  " Different numbers of quadrature points!");
271  homemade_assert_msg(lambda_corr_size == phi_corrected_dof,
272  " Incompatible corrected shape table and barycentric coordinates!");
273  homemade_assert_msg(lambda_inter_size == phi_inter_dof,
274  " Incompatible intersection shape table and barycentric coordinates!");
275 
276  // phi_A,i (qp) = l_i,1 * phi_I,1 (qp) + l_i,2 * phi_I,2 (qp) ...
277  // = sum [ l_i,j * phi_I,j (qp) ]
278  for(unsigned int qp = 0; qp < phi_inter_qp; ++qp)
279  {
280  for(unsigned int iii = 0; iii < phi_corrected_dof; ++iii)
281  {
282  dphi_corrected[iii][qp] = 0;
283  for(unsigned int jjj = 0; jjj < phi_inter_dof; ++jjj)
284  {
285  dphi_corrected[iii][qp] += lambda_weights[iii][jjj]*dphi_inter[jjj][qp];
286  }
287  }
288  }
289  };
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
void carl::assemble_coupling_matrices::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 at line 209 of file assemble_coupling.cpp.

212  {
213  unsigned int lambda_corr_size = lambda_weights.size();
214  unsigned int lambda_inter_size = lambda_weights[0].size();
215 
216  unsigned int phi_inter_dof = phi_inter.size();
217  unsigned int phi_inter_qp = phi_inter[0].size();
218 
219  unsigned int phi_corrected_dof = phi_corrected.size();
220  unsigned int phi_corrected_qp = phi_corrected[0].size();
221 
222  /*
223  * lambda : [ n_dofs_corr ] x [ n_dofs_inter ]
224  * inter : [ n_dofs_inter ] x [ n_qp ]
225  * corr : [ n_dofs_corr ] x [ n_qp ]
226  */
227 
228  homemade_assert_msg(phi_inter_qp == phi_corrected_qp,
229  " Different numbers of quadrature points!");
230  homemade_assert_msg(lambda_corr_size == phi_corrected_dof,
231  " Incompatible corrected shape table and barycentric coordinates!");
232  homemade_assert_msg(lambda_inter_size == phi_inter_dof,
233  " Incompatible intersection shape table and barycentric coordinates!");
234 
235  // phi_A,i (qp) = l_i,1 * phi_I,1 (qp) + l_i,2 * phi_I,2 (qp) ...
236  // = sum [ l_i,j * phi_I,j (qp) ]
237  for(unsigned int qp = 0; qp < phi_inter_qp; ++qp)
238  {
239  for(unsigned int iii = 0; iii < phi_corrected_dof; ++iii)
240  {
241  phi_corrected[iii][qp] = 0;
242  for(unsigned int jjj = 0; jjj < phi_inter_dof; ++jjj)
243  {
244  phi_corrected[iii][qp] += lambda_weights[iii][jjj]*phi_inter[jjj][qp];
245  }
246  }
247  }
248  }
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
void carl::assemble_coupling_matrices::set_coupling_parameters ( const std::string &  name,
double  coupling_rigidity,
double  coupling_width 
)

Set physical parameters.

Definition at line 202 of file assemble_coupling.cpp.

204  {
205  m_coupling_rigidityMap[name] = coupling_rigidity;
206  m_coupling_widthMap[name] = coupling_width;
207  };
carl::assemble_coupling_matrices::m_coupling_widthMap
std::map< std::string, double > m_coupling_widthMap
Definition: assemble_coupling.h:278
carl::assemble_coupling_matrices::m_coupling_rigidityMap
std::map< std::string, double > m_coupling_rigidityMap
Definition: assemble_coupling.h:277
libMesh::EquationSystems & carl::assemble_coupling_matrices::set_Restricted_BIG_EquationSystem ( const std::string &  name,
libMesh::MeshBase &  R_BIGMesh 
)

Definition at line 120 of file assemble_coupling.cpp.

121  {
122  m_bHasDefinedMeshRestrictions = true;
123  libMesh::EquationSystems* EqSystemPtr = new libMesh::EquationSystems(
124  R_BIGMesh);
125 
126  m_R_BIG_EquationSystem.first = name;
127  m_R_BIG_EquationSystem.second = EqSystemPtr;
128 
129  return *EqSystemPtr;
130  };
carl::assemble_coupling_matrices::m_bHasDefinedMeshRestrictions
bool m_bHasDefinedMeshRestrictions
Definition: assemble_coupling.h:274
carl::assemble_coupling_matrices::m_R_BIG_EquationSystem
std::pair< std::string, libMesh::EquationSystems * > m_R_BIG_EquationSystem
Definition: assemble_coupling.h:256
void carl::assemble_coupling_matrices::use_H1_coupling ( std::string  name)

Definition at line 406 of file assemble_coupling.cpp.

407  {
408  m_bUseH1Coupling[name] = true;
409  };
carl::assemble_coupling_matrices::m_bUseH1Coupling
std::map< std::string, bool > m_bUseH1Coupling
Definition: assemble_coupling.h:270
void carl::assemble_coupling_matrices::use_L2_coupling ( std::string  name)

Definition at line 411 of file assemble_coupling.cpp.

412  {
413  m_bUseH1Coupling[name] = false;
414  };
carl::assemble_coupling_matrices::m_bUseH1Coupling
std::map< std::string, bool > m_bUseH1Coupling
Definition: assemble_coupling.h:270

Member Data Documentation

bool carl::assemble_coupling_matrices::m_bHasAssembled_BIG
protected

Definition at line 272 of file assemble_coupling.h.

std::map<std::string, bool> carl::assemble_coupling_matrices::m_bHasAssembled_micro
protected

Definition at line 271 of file assemble_coupling.h.

bool carl::assemble_coupling_matrices::m_bHasDefinedMeshRestrictions
protected

Definition at line 274 of file assemble_coupling.h.

std::pair<std::string, libMesh::EquationSystems*> carl::assemble_coupling_matrices::m_BIG_EquationSystem
protected

Definition at line 255 of file assemble_coupling.h.

std::map<std::string, bool> carl::assemble_coupling_matrices::m_bUseH1Coupling
protected

Definition at line 270 of file assemble_coupling.h.

bool carl::assemble_coupling_matrices::m_bUseNullSpace_BIG
protected

Definition at line 273 of file assemble_coupling.h.

const libMesh::Parallel::Communicator* carl::assemble_coupling_matrices::m_comm
protected

Definition at line 252 of file assemble_coupling.h.

std::map<std::string, double> carl::assemble_coupling_matrices::m_coupling_rigidityMap
protected

Definition at line 277 of file assemble_coupling.h.

std::map<std::string, double> carl::assemble_coupling_matrices::m_coupling_widthMap
protected

Definition at line 278 of file assemble_coupling.h.

std::map<std::string, libMesh::PetscMatrix<libMesh::Number>*> carl::assemble_coupling_matrices::m_couplingMatrixMap_mediator_BIG
protected

Definition at line 266 of file assemble_coupling.h.

std::map<std::string, libMesh::PetscMatrix<libMesh::Number>*> carl::assemble_coupling_matrices::m_couplingMatrixMap_mediator_mediator
protected

Definition at line 267 of file assemble_coupling.h.

std::map<std::string, libMesh::PetscMatrix<libMesh::Number>*> carl::assemble_coupling_matrices::m_couplingMatrixMap_mediator_micro
protected

Definition at line 265 of file assemble_coupling.h.

std::map<std::string, libMesh::EquationSystems*> carl::assemble_coupling_matrices::m_inter_EquationSystemMap
protected

Definition at line 261 of file assemble_coupling.h.

std::map<std::string, libMesh::EquationSystems*> carl::assemble_coupling_matrices::m_mediator_EquationSystemMap
protected

Definition at line 262 of file assemble_coupling.h.

std::map<std::string, libMesh::EquationSystems*> carl::assemble_coupling_matrices::m_micro_EquationSystemMap
protected

Definition at line 258 of file assemble_coupling.h.

std::pair<std::string, libMesh::EquationSystems*> carl::assemble_coupling_matrices::m_R_BIG_EquationSystem
protected

Definition at line 256 of file assemble_coupling.h.

std::map<std::string, libMesh::EquationSystems*> carl::assemble_coupling_matrices::m_R_micro_EquationSystemMap
protected

Definition at line 259 of file assemble_coupling.h.

The documentation for this class was generated from the following files: