CArl
Code Arlequin / C++ implementation
Classes | Enumerations | Functions
carl Namespace Reference

Classes

class  assemble_coupling_matrices
 
struct  coupling_assemble_coupling_input_params
 Structure containing the parameters for the construction of the coupling matrices. More...
 
class  coupling_matrices_3
 
struct  feti_iterate_params
 Structure containing the parameters for the setup initialization of the FETI solver. More...
 
class  FETI_Operations
 Class containing the operations needed for the FETI solver. More...
 
struct  feti_set_sol_params
 Structure containing the parameters for the setup initialization of the FETI solver. More...
 
struct  feti_setup_finish_params
 Structure containing the parameters for the setup initialization of the FETI solver. More...
 
struct  feti_setup_init_params
 Structure containing the parameters for the setup initialization of the FETI solver. More...
 
class  Intersection_Search
 Class containing the structure needed to find all the intersections between two meshes, inside the coupling region mesh. More...
 
class  Intersection_Tools
 Class with a series of methods to find the intersections between libMesh's elements, using CGAL internally. More...
 
struct  IntersectionData
 
struct  libmesh_apply_solution_input_params
 
class  libMesh_fe_addresses_3
 
struct  libmesh_solve_linear_system_input_params
 
class  Mesh_Intersection
 Class used to construct the intersection meshes. More...
 
class  Mesh_restriction
 Class used to build a restriction of a parent mesh to the coupling region. More...
 
struct  parallel_intersection_params
 Structure containing the parameters for the parallel intersection search test program (source: CArl_build_intersections.cpp) More...
 
class  Patch_construction
 Class used to build a mesh patch from a parent mesh and an coupling mesh element. More...
 
struct  PointHash_3D
 
struct  PointHash_3D_Equal
 
class  Solver_Files_Setup
 
class  Stitch_Meshes
 Class used to stitch together different meshes. More...
 
class  weight_parameter_function
 

Enumerations

enum  ClusterSchedulerType { LOCAL = 0, PBS = 1, SLURM = 2 }
 
enum  ExtSolverType { LIBMESH_LINEAR = 0, DUMMY = 1 }
 
enum  MediatorType { USE_MACRO = 0, USE_MICRO = 1, USE_EXTERNAL = 2 }
 
enum  BaseCGPrecondType { NO_PRECONDITIONER = 0, COUPLING_OPERATOR = 1, COUPLING_JACOBI = 2 }
 
enum  IterationStatus { ITERATING = 0, DIVERGED = 1, CONVERGED = 2 }
 
enum  IntersectionMeshingMethod { LIBMESH_TETGEN = 0, CGAL = 1 }
 
enum  SearchMethod { BRUTE = 0, FRONT = 1, BOTH = 2 }
 
enum  RBModesSystem { MACRO = 0, MICRO = 1 }
 

Functions

void get_assemble_coupling_input_params (GetPot &field_parser, coupling_assemble_coupling_input_params &input_params)
 Parser function for the construction of the coupling matrices. More...
 
void get_input_params (GetPot &field_parser, feti_iterate_params &input_params)
 Parser function for the coupled solver test programs. More...
 
void get_input_params (GetPot &field_parser, feti_set_sol_params &input_params)
 Parser function for the coupled solver test programs. More...
 
void get_input_params (GetPot &field_parser, feti_setup_finish_params &input_params)
 Parser function for the coupled solver test programs. More...
 
void get_input_params (GetPot &field_parser, feti_setup_init_params &input_params)
 Parser function for the coupled solver test programs. More...
 
void get_intersection_input_params (GetPot &field_parser, parallel_intersection_params &input_params)
 Parser function for the parallel intersection search program (source: CArl_build_intersections.cpp) More...
 
void get_input_params (GetPot &field_parser, libmesh_solve_linear_system_input_params &input_params)
 Parser function for the coupled solver test programs. More...
 
void print_input_params (const std::string &output_filename, libmesh_solve_linear_system_input_params &input_params)
 Function used to generate a solver input file from "input_params". More...
 
void get_input_params (GetPot &field_parser, libmesh_apply_solution_input_params &input_params)
 Parser function for mesh deformation (apply solution) programs. More...
 
std::string ClusterSchedulerType_to_string (ClusterSchedulerType input)
 
std::string BaseCGPrecondType_to_string (BaseCGPrecondType input)
 
std::string ExtSolverType_to_string (ExtSolverType input)
 
std::string exec_command (const std::string &cmd)
 
void invert_index_unordered_map (const std::unordered_map< int, int > &input_map, std::unordered_map< int, int > &output_map)
 
template<typename T >
void jump_lines (T &filestream, unsigned int numberOfLines=1)
 
int voigt_index_converter (int aaa, int bbb)
 
void print_stats_to_file (std::vector< double > &vec_data, const std::string filename)
 
template<typename Sys >
void reduced_system_init (Sys &system_input)
 
void set_weight_function_domain_idx (std::string &filename, int &domain_Idx_BIG, int &nb_of_domain_Idx_micro, std::vector< int > &domain_Idx_micro, std::vector< int > &domain_Idx_coupling)
 
void build_intersection_and_restriction_tables (const libMesh::Parallel::Communicator &WorldComm, const std::string &intersection_full_table_Filename, const std::string &equivalence_table_A_Filename, const std::string &equivalence_table_B_Filename, std::vector< carl::IntersectionData > &intersection_full_table, std::unordered_map< int, int > &equivalence_table_A_to_R_A, std::unordered_map< int, int > &equivalence_table_B_to_R_B, std::unordered_map< int, int > &equivalence_table_R_A_to_A, std::unordered_map< int, int > &equivalence_table_R_B_to_B)
 
void generate_intersection_tables_partial (std::string &intersection_table_restrict_B_Filename, std::string &intersection_table_I_Filename, std::unordered_map< int, int > &mesh_restrict_ElemMap, std::unordered_map< int, int > &mesh_micro_ElemMap, std::unordered_map< int, int > &mesh_inter_ElemMap, std::vector< std::pair< int, int > > &intersection_table_restrict_B, std::unordered_multimap< int, int > &intersection_table_I)
 
void generate_intersection_tables_full (std::string &equivalence_table_restrict_A_Filename, std::string &intersection_table_restrict_B_Filename, std::string &intersection_table_I_Filename, std::unordered_map< int, int > &mesh_restrict_ElemMap, std::unordered_map< int, int > &mesh_micro_ElemMap, std::unordered_map< int, int > &mesh_BIG_ElemMap, std::unordered_map< int, int > &mesh_inter_ElemMap, std::unordered_map< int, int > &equivalence_table_restrict_A, std::vector< std::pair< int, int > > &intersection_table_restrict_B, std::unordered_multimap< int, int > &intersection_table_I)
 
void set_equivalence_tables (const libMesh::Parallel::Communicator &WorldComm, const std::string &equivalence_table_A_Filename, const std::string &equivalence_table_B_Filename, std::unordered_map< int, int > &equivalence_table_A_to_R_A, std::unordered_map< int, int > &equivalence_table_B_to_R_B, std::unordered_map< int, int > &equivalence_table_R_A_to_A, std::unordered_map< int, int > &equivalence_table_R_B_to_B)
 
void set_restricted_intersection_pairs_table (const std::unordered_map< int, std::pair< int, int > > &full_intersection_pairs_map, const std::unordered_map< int, int > &equivalence_table_A_to_R_A, const std::unordered_map< int, int > &equivalence_table_B_to_R_B, std::unordered_map< int, std::pair< int, int > > &full_intersection_restricted_pairs_map)
 
void set_full_intersection_tables (const libMesh::Parallel::Communicator &WorldComm, const std::string &intersection_full_table_Filename, std::unordered_map< int, std::pair< int, int > > &full_intersection_pairs_map, std::unordered_map< int, int > &full_intersection_meshI_to_inter_map)
 
void set_intersection_tables (const libMesh::Parallel::Communicator &WorldComm, const libMesh::Mesh &mesh_intersection, const std::string &intersection_full_table_Filename, const std::string &equivalence_table_A_Filename, const std::string &equivalence_table_B_Filename, const std::unordered_map< int, int > &equivalence_table_A_to_R_A, const std::unordered_map< int, int > &equivalence_table_B_to_R_B, std::unordered_map< int, std::pair< int, int > > &full_intersection_pairs_map, std::unordered_map< int, std::pair< int, int > > &full_intersection_restricted_pairs_map, std::unordered_map< int, int > &local_intersection_meshI_to_inter_map)
 
void read_local_intersection_tables (const libMesh::Parallel::Communicator &WorldComm, const std::string &intersection_local_table_Filename, std::unordered_map< int, std::pair< int, int > > &local_intersection_pairs_map, std::unordered_map< int, int > &local_intersection_meshI_to_inter_map)
 
void set_local_intersection_tables (const libMesh::Parallel::Communicator &WorldComm, const libMesh::Mesh &mesh_intersection, const std::string &intersection_local_table_Filename, const std::string &equivalence_table_A_Filename, const std::string &equivalence_table_B_Filename, const std::unordered_map< int, int > &equivalence_table_A_to_R_A, const std::unordered_map< int, int > &equivalence_table_B_to_R_B, std::unordered_map< int, std::pair< int, int > > &local_intersection_pairs_map, std::unordered_map< int, std::pair< int, int > > &local_intersection_restricted_pairs_map, std::unordered_map< int, int > &local_intersection_meshI_to_inter_map)
 
void set_global_mediator_system_intersection_lists (const libMesh::Parallel::Communicator &WorldComm, const std::string &intersection_global_table_Filename, const std::unordered_map< int, int > &equivalence_table_system_to_mediator, const std::unordered_map< int, int > &equivalence_table_mediator_to_system, std::unordered_multimap< int, int > &inter_mediator_A, std::unordered_multimap< int, int > &inter_mediator_B)
 
void repartition_system_meshes (const libMesh::Parallel::Communicator &WorldComm, libMesh::Mesh &mesh_input, libMesh::Mesh &mesh_intersect, std::unordered_map< int, std::pair< int, int > > &local_intersection_pairs_map, bool bUseSecond=true)
 
void broadcast_index_unordered_map (std::unordered_map< int, int > &index_map, const libMesh::Parallel::Communicator &CommComm, int origin_rank=0)
 
template<typename T >
void MPI_reduce_vector (std::vector< T > &r, int root, const libMesh::Parallel::Communicator &Comm)
 
void lump_matrix (libMesh::PetscMatrix< libMesh::Number > &matrixInput, libMesh::PetscMatrix< libMesh::Number > &matrixOutput)
 
void lump_matrix_and_invert (libMesh::PetscMatrix< libMesh::Number > &matrixInput, libMesh::PetscMatrix< libMesh::Number > &matrixOutput)
 
void lump_matrix_and_invert (libMesh::PetscMatrix< libMesh::Number > &matrixInput, libMesh::PetscVector< libMesh::Number > &vecOutput)
 
void print_matrix (libMesh::PetscMatrix< libMesh::Number > &CouplingTestMatrix)
 
void print_matrix_dim (libMesh::PetscMatrix< libMesh::Number > &CouplingTestMatrix, bool bDetailed=false)
 
void print_matrix_info (libMesh::PetscMatrix< libMesh::Number > &InputMatrix, std::ostream &os=libMesh::out)
 
void print_matrix_col_line_sum (libMesh::PetscMatrix< libMesh::Number > &CouplingTestMatrix, const std::string name_base)
 
void print_matrix_matlab (libMesh::PetscMatrix< libMesh::Number > &CouplingTestMatrix, const std::string name_base)
 
void solve_linear_PETSC (libMesh::PetscMatrix< libMesh::Number > &A, libMesh::PetscVector< libMesh::Number > &b, libMesh::PetscVector< libMesh::Number > &x, KSP &ksp, PC &pc)
 
void check_coupling_matrix (libMesh::PetscMatrix< libMesh::Number > &CouplingTestMatrix, libMesh::Mesh &IntersectionMesh, libMesh::Real CouplingScale, const std::string matrixtype, int n_var=3)
 
void write_PETSC_vector (libMesh::PetscVector< libMesh::Number > &input_vec, const std::string &filename, int dim=1)
 
void read_PETSC_vector (libMesh::PetscVector< libMesh::Number > &input_vec, const std::string &filename)
 
void print_PETSC_vector (libMesh::PetscVector< libMesh::Number > &input_vec, const std::string &filename)
 
void write_PETSC_vector (Vec input_vec, const std::string &filename, int rank, MPI_Comm comm=PETSC_COMM_WORLD, int dim=1)
 
void read_PETSC_vector (Vec input_vec, const std::string &filename, MPI_Comm comm=PETSC_COMM_WORLD)
 
void write_PETSC_matrix (Mat input_mat, const std::string &filename, int rank, MPI_Comm comm=PETSC_COMM_WORLD, int dim=1)
 
void write_PETSC_matrix (libMesh::PetscMatrix< libMesh::Number > &input_mat, const std::string &filename, int dim=1)
 
void read_PETSC_matrix (Mat input_mat, const std::string &filename, MPI_Comm comm=PETSC_COMM_WORLD)
 
void read_PETSC_matrix (libMesh::PetscMatrix< libMesh::Number > &input_mat, const std::string &filename)
 
void write_PETSC_vector_MATLAB (Vec input_vec, const std::string &filename, MPI_Comm comm=PETSC_COMM_WORLD)
 
void write_PETSC_matrix_MATLAB (Mat input_mat, const std::string &filename, MPI_Comm comm=PETSC_COMM_WORLD)
 
void attach_rigid_body_mode_vectors (libMesh::PetscMatrix< libMesh::Number > &mat_sys, const std::string &filename_base, int nb_of_vecs, int dimension)
 
void create_PETSC_dense_matrix_from_vectors (const Vec *vecs_in, int nb_vecs, Mat &matrix_out)
 
void PETSC_invert_dense_matrix (Mat &matrix_in, Mat &matrix_out)
 
void PETSC_MatMultScale_Bcast (Mat mat_seq, Vec vec_seq_in, Vec vec_seq_out, double a_const)
 

Enumeration Type Documentation

enum carl::BaseCGPrecondType
Enumerator
NO_PRECONDITIONER 
COUPLING_OPERATOR 
COUPLING_JACOBI 

Definition at line 35 of file common_enums.h.

35  {
36  NO_PRECONDITIONER = 0, // Identity matrix
37  COUPLING_OPERATOR = 1, // C_RR
38  COUPLING_JACOBI = 2 // diagonal(C_RR)
39 };
enum carl::ClusterSchedulerType
Enumerator
LOCAL 
PBS 
SLURM 

Definition at line 14 of file common_enums.h.

14  {
15  LOCAL = 0, // No scheduler present, will use carl::exec_command
16  PBS = 1, // PBS / Torque
17  SLURM = 2 // SLURM, not implemented right now
18 };
enum carl::ExtSolverType
Enumerator
LIBMESH_LINEAR 
DUMMY 

Definition at line 22 of file common_enums.h.

22  {
23  LIBMESH_LINEAR = 0,
24  DUMMY = 1
25 };
enum carl::IntersectionMeshingMethod
Enumerator
LIBMESH_TETGEN 
CGAL 

Definition at line 47 of file common_enums.h.

47  {
48  LIBMESH_TETGEN = 0, // libMesh Tetgen algorithm, problematic with
49  // Intel compilers
50  CGAL = 1 // Intersection meshing algorithm using
51  // CGAL's Triangulation_3
52 };
enum carl::IterationStatus
Enumerator
ITERATING 
DIVERGED 
CONVERGED 

Definition at line 41 of file common_enums.h.

41  {
42  ITERATING = 0,
43  DIVERGED = 1,
44  CONVERGED = 2
45 };
enum carl::MediatorType
Enumerator
USE_MACRO 
USE_MICRO 
USE_EXTERNAL 

Definition at line 28 of file common_enums.h.

28  {
29  USE_MACRO = 0,
30  USE_MICRO = 1,
31  USE_EXTERNAL = 2
32 };
enum carl::RBModesSystem
Enumerator
MACRO 
MICRO 

Definition at line 61 of file common_enums.h.

61  {
62  MACRO = 0,
63  MICRO = 1
64 };
MACRO
Definition: ext_solver_libmesh_enums.h:16
MICRO
Definition: ext_solver_libmesh_enums.h:17
enum carl::SearchMethod
Enumerator
BRUTE 
FRONT 
BOTH 

Definition at line 54 of file common_enums.h.

55 {
56  BRUTE = 0,
57  FRONT = 1,
58  BOTH = 2
59 };

Function Documentation

void carl::attach_rigid_body_mode_vectors ( libMesh::PetscMatrix< libMesh::Number > &  mat_sys,
const std::string &  filename_base,
int  nb_of_vecs,
int  dimension 
)

Definition at line 383 of file PETSC_matrix_operations.cpp.

385 {
386  // Read the dummy vector
387  Vec rb_vecs[6];
388 
389  PetscInt local_N;
390  MatGetLocalSize(mat_sys.mat(),NULL,&local_N);
391 
392  std::string vector_path;
393 
394  // Read the vectors
395  for(int iii = 0; iii < nb_of_vecs; ++iii)
396  {
397  VecCreate(mat_sys.comm().get(),&rb_vecs[iii]);
398  VecSetSizes(rb_vecs[iii],local_N,mat_sys.n());
399 
400  vector_path = filename_base + "_" + std::to_string(iii) + "_n_" + std::to_string(nb_of_vecs) + ".petscvec";
401  read_PETSC_vector(rb_vecs[iii],vector_path,mat_sys.comm().get());
402  }
403 
404  MatNullSpace nullsp_sys;
405  MatNullSpaceCreate(mat_sys.comm().get(), PETSC_FALSE, nb_of_vecs, rb_vecs, &nullsp_sys);
406  MatSetNullSpace(mat_sys.mat(),nullsp_sys);
407 
408  MatNullSpaceDestroy(&nullsp_sys);
409  for(int iii = 0; iii < nb_of_vecs; ++iii)
410  {
411  VecDestroy(&rb_vecs[iii]);
412  }
413 };
carl::read_PETSC_vector
void read_PETSC_vector(libMesh::PetscVector< libMesh::Number > &input_vec, const std::string &filename)
Definition: PETSC_matrix_operations.cpp:377
std::string carl::BaseCGPrecondType_to_string ( BaseCGPrecondType  input)

Definition at line 27 of file common_functions.cpp.

28 {
29  switch (input)
30  {
31  case BaseCGPrecondType::NO_PRECONDITIONER : return "NONE";
32  break;
33 
34  case BaseCGPrecondType::COUPLING_OPERATOR : return "Coupling_operator";
35  break;
36 
37  case BaseCGPrecondType::COUPLING_JACOBI : return "Coupling_operator_jacobi";
38  break;
39  }
40 
41  homemade_error_msg("Invalid enumerate argument!");
42  return "";
43 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::broadcast_index_unordered_map ( std::unordered_map< int, int > &  index_map,
const libMesh::Parallel::Communicator &  CommComm,
int  origin_rank = 0 
)

Definition at line 3 of file mpi_carl_tools.cpp.

7 {
8  int rank = CommComm.rank();
9  int dummy_vector_size = -1;
10  std::vector<int> dummy_vector;
11 
12  if(rank == origin_rank)
13  {
14  dummy_vector_size = index_map.size();
15  dummy_vector.resize(2*dummy_vector_size);
16 
17  std::unordered_map<int,int>::iterator mapIt = index_map.begin();
18  std::unordered_map<int,int>::iterator end_mapIt = index_map.end();
19  int dummy_iii = 0;
20  for( ; mapIt != end_mapIt; ++mapIt)
21  {
22  dummy_vector[2*dummy_iii] = mapIt->first;
23  dummy_vector[2*dummy_iii + 1] = mapIt->second;
24  ++dummy_iii;
25  }
26  }
27 
28  CommComm.barrier();
29  CommComm.broadcast(dummy_vector_size,origin_rank);
30 
31  if( rank != origin_rank )
32  {
33  dummy_vector.resize(2*dummy_vector_size);
34  }
35  CommComm.broadcast(dummy_vector,origin_rank);
36 
37  if( rank != origin_rank)
38  {
39  index_map.reserve(dummy_vector_size);
40  for(int iii = 0; iii < dummy_vector_size; ++iii)
41  {
42  index_map[dummy_vector[2*iii]] = dummy_vector[2*iii+1];
43  }
44  }
45 };
void carl::build_intersection_and_restriction_tables ( const libMesh::Parallel::Communicator &  WorldComm,
const std::string &  intersection_full_table_Filename,
const std::string &  equivalence_table_A_Filename,
const std::string &  equivalence_table_B_Filename,
std::vector< carl::IntersectionData > &  intersection_full_table,
std::unordered_map< int, int > &  equivalence_table_A_to_R_A,
std::unordered_map< int, int > &  equivalence_table_B_to_R_B,
std::unordered_map< int, int > &  equivalence_table_R_A_to_A,
std::unordered_map< int, int > &  equivalence_table_R_B_to_B 
)

Definition at line 182 of file mesh_tables.cpp.

193 {
194  int rank = WorldComm.rank();
195 
196  // While the equivalence tables are saved as unordered maps, it's easier to
197  // save them as vectors at first, broadcast them, and them reconvert to maps
198  std::vector<int> dummy_equivalence_table_A;
199  std::vector<int> dummy_equivalence_table_B;
200  std::vector<int> dummy_intersection_full_table;
201 
202  int nbOfIntersections = -1;
203  int nbOfInterElems = -1;
204  int nbOfRestricted_A_Elems = -1;
205  int nbOfRestricted_B_Elems = -1;
206 
207  // Do the file reading job on proc 0
208  if(rank == 0)
209  {
210  std::ifstream intersection_full_file(intersection_full_table_Filename);
211 
212  intersection_full_file >> nbOfIntersections >> nbOfInterElems;
213  dummy_intersection_full_table.resize(4*nbOfInterElems);
214 
215  for(int iii = 0; iii < nbOfInterElems; ++iii)
216  {
217  intersection_full_file
218  >> dummy_intersection_full_table[4*iii]
219  >> dummy_intersection_full_table[4*iii + 1]
220  >> dummy_intersection_full_table[4*iii + 2]
221  >> dummy_intersection_full_table[4*iii + 3];
222  }
223  intersection_full_file.close();
224 
225  std::ifstream equivalence_A_file(equivalence_table_A_Filename);
226 
227  equivalence_A_file >> nbOfRestricted_A_Elems;
228  dummy_equivalence_table_A.resize(2*nbOfRestricted_A_Elems);
229 
230  for(int iii = 0; iii < nbOfRestricted_A_Elems; ++iii)
231  {
232  equivalence_A_file >> dummy_equivalence_table_A[2*iii]
233  >> dummy_equivalence_table_A[2*iii + 1];
234  }
235  equivalence_A_file.close();
236 
237  std::ifstream equivalence_B_file(equivalence_table_B_Filename);
238 
239  equivalence_B_file >> nbOfRestricted_B_Elems;
240  dummy_equivalence_table_B.resize(2*nbOfRestricted_B_Elems);
241 
242  for(int iii = 0; iii < nbOfRestricted_B_Elems; ++iii)
243  {
244  equivalence_B_file >> dummy_equivalence_table_B[2*iii]
245  >> dummy_equivalence_table_B[2*iii + 1];
246  }
247  equivalence_B_file.close();
248  }
249 
250  // Broadcast the sizes and resize on other procs
251  WorldComm.broadcast(nbOfIntersections);
252  WorldComm.broadcast(nbOfInterElems);
253  WorldComm.broadcast(nbOfRestricted_A_Elems);
254  WorldComm.broadcast(nbOfRestricted_B_Elems);
255 
256  if(rank != 0)
257  {
258  dummy_intersection_full_table.resize(4*nbOfInterElems);
259  dummy_equivalence_table_A.resize(2*nbOfRestricted_A_Elems);
260  dummy_equivalence_table_B.resize(2*nbOfRestricted_B_Elems);
261  }
262 
263  WorldComm.barrier();
264 
265  WorldComm.broadcast(dummy_intersection_full_table);
266  WorldComm.broadcast(dummy_equivalence_table_A);
267  WorldComm.broadcast(dummy_equivalence_table_B);
268 
269  // Convert back to unoredered maps and pair vectors
270  intersection_full_table.resize(nbOfInterElems);
271  equivalence_table_A_to_R_A.reserve(nbOfRestricted_A_Elems);
272  equivalence_table_B_to_R_B.reserve(nbOfRestricted_B_Elems);
273 
274  equivalence_table_R_A_to_A.reserve(nbOfRestricted_A_Elems);
275  equivalence_table_R_B_to_B.reserve(nbOfRestricted_B_Elems);
276 
277  for(int iii = 0; iii < nbOfInterElems; ++iii)
278  {
279  intersection_full_table[iii].InterMeshIdx
280  = dummy_intersection_full_table[4*iii];
281  intersection_full_table[iii].AMeshIdx
282  = dummy_intersection_full_table[4*iii + 1];
283  intersection_full_table[iii].BMeshIdx
284  = dummy_intersection_full_table[4*iii + 2];
285  intersection_full_table[iii].IntersectionID
286  = dummy_intersection_full_table[4*iii + 3];
287  }
288 
289  for(int iii = 0; iii < nbOfRestricted_A_Elems; ++iii)
290  {
291  equivalence_table_R_A_to_A[dummy_equivalence_table_A[2*iii]] =
292  dummy_equivalence_table_A[2*iii + 1];
293  equivalence_table_A_to_R_A[dummy_equivalence_table_A[2*iii + 1]] =
294  dummy_equivalence_table_A[2*iii];
295  }
296 
297  for(int iii = 0; iii < nbOfRestricted_B_Elems; ++iii)
298  {
299  equivalence_table_R_B_to_B[dummy_equivalence_table_B[2*iii]] =
300  dummy_equivalence_table_B[2*iii + 1];
301  equivalence_table_B_to_R_B[dummy_equivalence_table_B[2*iii + 1]] =
302  dummy_equivalence_table_B[2*iii];
303  }
304 };
void carl::check_coupling_matrix ( libMesh::PetscMatrix< libMesh::Number > &  CouplingTestMatrix,
libMesh::Mesh &  IntersectionMesh,
libMesh::Real  CouplingScale,
const std::string  matrixtype,
int  n_var = 3 
)

Definition at line 284 of file PETSC_matrix_operations.cpp.

289 {
290  std::cout << "| " << matrixType << std::endl;
291  libMesh::Real accumulator = 0;
292 
293  std::cout << "| M_i,j : " << CouplingTestMatrix.m() << " x " << CouplingTestMatrix.n() << std::endl;
294  for(unsigned int iii = 0; iii < CouplingTestMatrix.m(); ++iii)
295  {
296  for(unsigned int jjj = 0; jjj < CouplingTestMatrix.n(); ++jjj)
297  {
298  accumulator += CouplingTestMatrix(iii,jjj);
299  }
300  }
301 
302  libMesh::Real vol = 0;
303  libMesh::Elem* silly_elem;
304  for(libMesh::MeshBase::element_iterator itBegin = IntersectionMesh.elements_begin();
305  itBegin != IntersectionMesh.elements_end();
306  ++itBegin)
307  {
308  silly_elem = *itBegin;
309  vol += silly_elem->volume();
310  }
311 
312  libMesh::Real difference = accumulator - n_var*CouplingScale*vol;
313 
314  std::cout << "|" << std::endl;
315  std::cout << "| Sum( M_i,j ) = " << accumulator << std::endl;
316  std::cout << "| n * C * Volume = " << n_var*CouplingScale*vol << std::endl;
317  std::cout << "| > Difference = " << difference << std::endl << std::endl;
318 }
std::string carl::ClusterSchedulerType_to_string ( ClusterSchedulerType  input)

Definition at line 9 of file common_functions.cpp.

10 {
11  switch (input)
12  {
13  case ClusterSchedulerType::LOCAL : return "LOCAL";
14  break;
15 
16  case ClusterSchedulerType::PBS : return "PBS";
17  break;
18 
19  case ClusterSchedulerType::SLURM : return "SLURM";
20  break;
21  }
22 
23  homemade_error_msg("Invalid enumerate argument!");
24  return "";
25 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::create_PETSC_dense_matrix_from_vectors ( const Vec *  vecs_in,
int  nb_vecs,
Mat &  matrix_out 
)

Definition at line 456 of file PETSC_matrix_operations.cpp.

457 {
458  // Get the vectors' dimensions and create the matrix
459  PetscInt M_local, M;
460  VecGetLocalSize(vecs_in[0],&M_local);
461  VecGetSize(vecs_in[0],&M);
462 
463  MatCreateDense(PETSC_COMM_WORLD,M_local,PETSC_DECIDE,M,nb_vecs,NULL,&matrix_out);
464 
465  // Get the ownership ranges and the row indexes
466  PetscInt own_low, own_high;
467  VecGetOwnershipRange(vecs_in[0],&own_low,&own_high);
468 
469  std::vector<PetscInt> row(M_local,0);
470 
471  for(int iii = 0; iii < M_local; ++iii)
472  {
473  row[iii] = own_low + iii;
474  }
475 
476  // Set the values
477  PetscScalar* data;
478  for(int iii = 0; iii < nb_vecs; ++iii)
479  {
480  VecGetArray(vecs_in[iii],&data);
481  MatSetValues(matrix_out,M_local,row.data(),1,&iii,data,INSERT_VALUES);
482  VecRestoreArray(vecs_in[iii],&data);
483  }
484 
485  // Assembly
486  MatAssemblyBegin(matrix_out,MAT_FINAL_ASSEMBLY);
487  MatAssemblyEnd(matrix_out,MAT_FINAL_ASSEMBLY);
488 }
std::string carl::exec_command ( const std::string &  cmd)

Definition at line 60 of file common_functions.cpp.

60  {
61  std::array<char, 128> buffer;
62  std::string result;
63  std::shared_ptr<FILE> pipe(popen(cmd.c_str(), "r"), pclose);
64  if (!pipe) throw std::runtime_error("popen() failed!");
65  while (!feof(pipe.get())) {
66  if (fgets(buffer.data(), 128, pipe.get()) != NULL)
67  result += buffer.data();
68  }
69  return result;
70 }
std::string carl::ExtSolverType_to_string ( ExtSolverType  input)

Definition at line 45 of file common_functions.cpp.

46 {
47  switch (input)
48  {
49  case ExtSolverType::LIBMESH_LINEAR : return "LIBMESH_LINEAR";
50  break;
51 
52  case ExtSolverType::DUMMY : return "DUMMY";
53  break;
54  }
55 
56  homemade_error_msg("Invalid enumerate argument!");
57  return "";
58 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::generate_intersection_tables_full ( std::string &  equivalence_table_restrict_A_Filename,
std::string &  intersection_table_restrict_B_Filename,
std::string &  intersection_table_I_Filename,
std::unordered_map< int, int > &  mesh_restrict_ElemMap,
std::unordered_map< int, int > &  mesh_micro_ElemMap,
std::unordered_map< int, int > &  mesh_BIG_ElemMap,
std::unordered_map< int, int > &  mesh_inter_ElemMap,
std::unordered_map< int, int > &  equivalence_table_restrict_A,
std::vector< std::pair< int, int > > &  intersection_table_restrict_B,
std::unordered_multimap< int, int > &  intersection_table_I 
)

Definition at line 104 of file mesh_tables.cpp.

115 {
116  std::ifstream table_restrict_A_file(equivalence_table_restrict_A_Filename);
117  std::ifstream table_restrict_B_file(intersection_table_restrict_B_Filename);
118  std::ifstream table_I_file(intersection_table_I_Filename);
119 
120  int nbOfEquivalences_restrict_A = -1;
121  int nbOfIntersections_restrict_B = -1;
122  int nbOfIntersectionsI = -1;
123  int nbOfTotalTetrasI = -1;
124  int dummyInt = -1;
125  int nbOfTetras = -1;
126  int tetraIdx = -1;
127 
128  int extIdxA = -1;
129  int extIdxB = -1;
130  int extIdxI = -1;
131  int extIdxR = -1;
132 
133  int idxRestrict = -1;
134  int idxA = -1;
135 
136  table_restrict_A_file >> nbOfEquivalences_restrict_A;
137  std::unordered_map<int,int> temp_equivalence_table_A_restrict(nbOfEquivalences_restrict_A);
138  equivalence_table_restrict_A.reserve(nbOfEquivalences_restrict_A);
139 
140  for(int iii = 0; iii < nbOfEquivalences_restrict_A; ++iii)
141  {
142  table_restrict_A_file >> extIdxR >> extIdxA;
143 
144  idxA = mesh_BIG_ElemMap[extIdxA];
145  idxRestrict = mesh_restrict_ElemMap[extIdxR];
146 
147  temp_equivalence_table_A_restrict[idxA] = idxRestrict;
148  equivalence_table_restrict_A[idxRestrict] = idxA;
149  }
150 
151  table_restrict_B_file >> nbOfIntersections_restrict_B;
152  table_I_file >> nbOfIntersectionsI >> nbOfTotalTetrasI;
153 
154  homemade_assert_msg(nbOfIntersections_restrict_B == nbOfIntersectionsI, "Incompatible intersection table files!");
155 
156  intersection_table_restrict_B.resize(nbOfIntersections_restrict_B);
157  intersection_table_I.reserve(nbOfTotalTetrasI);
158 
159  for(int iii = 0; iii < nbOfIntersections_restrict_B; ++iii)
160  {
161  table_restrict_B_file >> dummyInt >> extIdxA >> extIdxB;
162 
163  idxA = mesh_BIG_ElemMap[extIdxA];
164  intersection_table_restrict_B[iii].first = temp_equivalence_table_A_restrict[idxA];
165  intersection_table_restrict_B[iii].second = mesh_micro_ElemMap[extIdxB];
166 
167  table_I_file >> dummyInt >> nbOfTetras;
168  for(int jjj = 0; jjj < nbOfTetras; ++jjj)
169  {
170  table_I_file >> extIdxI;
171  tetraIdx = mesh_inter_ElemMap[extIdxI];
172  intersection_table_I.emplace(dummyInt,tetraIdx);
173  }
174  }
175 
176  table_restrict_A_file.close();
177  table_restrict_B_file.close();
178  table_I_file.close();
179 };
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
void carl::generate_intersection_tables_partial ( std::string &  intersection_table_restrict_B_Filename,
std::string &  intersection_table_I_Filename,
std::unordered_map< int, int > &  mesh_restrict_ElemMap,
std::unordered_map< int, int > &  mesh_micro_ElemMap,
std::unordered_map< int, int > &  mesh_inter_ElemMap,
std::vector< std::pair< int, int > > &  intersection_table_restrict_B,
std::unordered_multimap< int, int > &  intersection_table_I 
)

Definition at line 54 of file mesh_tables.cpp.

62 {
63  std::ifstream table_restrict_B_file(intersection_table_restrict_B_Filename);
64  std::ifstream table_I_file(intersection_table_I_Filename);
65 
66  int nbOfIntersections_restrict_B = -1;
67  int nbOfIntersectionsI = -1;
68  int nbOfTotalTetrasI = -1;
69  int dummyInt = -1;
70  int nbOfTetras = -1;
71  int tetraIdx = -1;
72 
73  int extIdxA = -1;
74  int extIdxB = -1;
75  int extIdxI = -1;
76 
77  table_restrict_B_file >> nbOfIntersections_restrict_B;
78  table_I_file >> nbOfIntersectionsI >> nbOfTotalTetrasI;
79 
80  homemade_assert_msg(nbOfIntersections_restrict_B == nbOfIntersectionsI, "Incompatible intersection table files!");
81 
82  intersection_table_restrict_B.resize(nbOfIntersections_restrict_B);
83  intersection_table_I.reserve(2*nbOfTotalTetrasI);
84 
85  for(int iii = 0; iii < nbOfIntersections_restrict_B; ++iii)
86  {
87  table_restrict_B_file >> dummyInt >> extIdxA >> extIdxB;
88  intersection_table_restrict_B[iii].first = mesh_restrict_ElemMap[extIdxA];
89  intersection_table_restrict_B[iii].second = mesh_micro_ElemMap[extIdxB];
90 
91  table_I_file >> dummyInt >> nbOfTetras;
92  for(int jjj = 0; jjj < nbOfTetras; ++jjj)
93  {
94  table_I_file >> extIdxI;
95  tetraIdx = mesh_inter_ElemMap[extIdxI];
96  intersection_table_I.emplace(dummyInt,tetraIdx);
97  }
98  }
99 
100  table_restrict_B_file.close();
101  table_I_file.close();
102 };
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
void carl::get_assemble_coupling_input_params ( GetPot &  field_parser,
coupling_assemble_coupling_input_params &  input_params 
)

Parser function for the construction of the coupling matrices.

Required parameters:

  • System and intersection meshes:
    • MeshA, -mA or --meshA : path to the mesh A.
    • MeshB, -mB or --meshB : path to the mesh B.
    • InterBase, -mI or --meshI : common path to the intersection meshes and tables.
  • Coupling parameters:
    • CouplingWidth or --ce : width of the coupling region (same unit as the meshes, $e$ in the $L_2$ coupling term).
    • CouplingRigidity or --ck : rigidity used for the coupling matrix (in MPa, if mm was used for the meshes, $\kappa$ in both the $L_2$ and $H_1$ terms).

Optional parameters:

  • Output:
    • OutputFolder or --output : base of the output folder. Default: "".
  • Restriction meshes and tables:
    • Mesh_A_Restriction, -mAR or --meshAR : path to the restricted mesh A (formed by elements of the mesh A intersecting the coupling region). Default: [InterBase]_A_restriction.msh.
    • Mesh_B_Restriction, -mBR or --meshBR : path to the restricted mesh B (formed by elements of the mesh A intersecting the coupling region). Default: [InterBase]_B_restriction.msh.
    • Mesh_A_RestrictionEquivalenceTable or --tableRA : path to the equivalence table between the mesh A and its restriction. Default: [InterBase]_A_restriction_restrict.dat.
    • Mesh_B_RestrictionEquivalenceTable or --tableRB : path to the equivalence table between the mesh B and its restriction. Default: [InterBase]_B_restriction_restrict.dat.
  • Mediator mesh:
    • MediatorMesh : choice of the mediator mesh. Values: UseRestricted_A or UseRestricted_B. Default: UseRestricted_A.

Definition at line 13 of file carl_assemble_coupling_input_parser.cpp.

14  {
15 
16  // Set mesh files
17  if (field_parser.search(3, "--meshA", "-mA", "MeshA")) {
18  input_params.mesh_BIG_file = field_parser.next(
19  input_params.mesh_BIG_file);
20  } else {
21  homemade_error_msg("Missing the A mesh file!");
22  }
23 
24  if (field_parser.search(3, "--meshB", "-mB", "MeshB")) {
25  input_params.mesh_micro_file = field_parser.next(
26  input_params.mesh_micro_file);
27  } else {
28  homemade_error_msg("Missing the B mesh file!");
29  }
30 
31  if (field_parser.search(3, "--meshI", "-mI", "InterBase")) {
32  input_params.common_inter_file = field_parser.next(
33  input_params.common_inter_file);
34  } else {
35  homemade_error_msg("Missing the path to the intersection files!");
36  }
37 
38  // Set coupling parameters
39  if( field_parser.search(2, "--ce","CouplingWidth") )
40  {
41  input_params.coupling_width = field_parser.next(input_params.coupling_width);
42  } else {
43  homemade_error_msg("Missing the coupling region width!");
44  }
45 
46  if( field_parser.search(2, "--ck","CouplingRigidity") )
47  {
48  input_params.coupling_rigidity = field_parser.next(input_params.coupling_rigidity);
49  } else {
50  homemade_error_msg("Missing the coupling rigidity!");
51  }
52 
53  // Output
54  if (field_parser.search(3, "--output", "-mO", "OutputFolder"))
55  {
56  input_params.output_folder = field_parser.next(
57  input_params.output_folder);
58  } else {
59  input_params.output_folder = "";
60  }
61 
62  if (field_parser.search(3, "--meshAR", "-mAR", "Mesh_A_Restriction")) {
63  input_params.mesh_restrict_BIG_file = field_parser.next(
64  input_params.mesh_restrict_BIG_file);
65  } else {
66  input_params.mesh_restrict_BIG_file = input_params.common_inter_file + "_A_restriction.msh";
67  }
68 
69  if (field_parser.search(3, "--meshBR", "-mBR", "Mesh_B_Restriction")) {
70  input_params.mesh_restrict_micro_file = field_parser.next(
71  input_params.mesh_restrict_micro_file);
72  } else {
73  input_params.mesh_restrict_micro_file = input_params.common_inter_file + "_B_restriction.msh";
74  }
75 
76  // Set the equivalence tables
77  if (field_parser.search(2, "--tableRA", "Mesh_A_RestrictionEquivalenceTable")) {
78  input_params.equivalence_table_restrict_BIG_file = field_parser.next(
79  input_params.equivalence_table_restrict_BIG_file);
80  } else {
81  input_params.equivalence_table_restrict_BIG_file = input_params.common_inter_file + "_A_restriction_restrict.dat";
82  }
83 
84  if (field_parser.search(2, "--tableRB", "Mesh_B_RestrictionEquivalenceTable")) {
85  input_params.equivalence_table_restrict_micro_file = field_parser.next(
86  input_params.equivalence_table_restrict_micro_file);
87  } else {
88  input_params.equivalence_table_restrict_micro_file = input_params.common_inter_file + "_B_restriction_restrict.dat";
89  }
90 
91  // Set the mediator mesh
92  input_params.mediator_type = carl::MediatorType::USE_MACRO;
93  input_params.mesh_mediator_file = input_params.mesh_restrict_BIG_file;
94 
95  std::string mediator_type;
96  if (field_parser.search(1,"MediatorMesh"))
97  {
98  mediator_type = field_parser.next(mediator_type);
99 
100  if(mediator_type == "UseRestricted_A")
101  {
102  input_params.mediator_type = carl::MediatorType::USE_MACRO;
103  input_params.mesh_mediator_file = input_params.mesh_restrict_BIG_file;
104  }
105  else if(mediator_type == "UseRestricted_B")
106  {
107  input_params.mediator_type = carl::MediatorType::USE_MICRO;
108  input_params.mesh_mediator_file = input_params.mesh_restrict_micro_file;
109  }
110  }
111 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::get_input_params ( GetPot &  field_parser,
libmesh_solve_linear_system_input_params &  input_params 
)

Parser function for the coupled solver test programs.

Required parameters:

  • SysMatrix : path to the system matrix file.
  • SysRHSVector : path to the system RHS vector file.
  • OutputBase : output filename base.

Optional parameters:

  • SysEps : relative convergence parameter. Default: 1e-5.
  • SysIterDiv : maximum number of iterations. Default: 1e3.
  • RBVectorBase : filename base to the rigid body mode vectors. The program expects that these vectors will be named as [RBVectorBase]_rb_vector_XYZ_n_[NbOfRBVectors].petscvec, where XYZ is an integer index going from 0 to NbOfRBVectors - 1. Setting this parameter sets bUseRBVectors to true - else, it is set to false.
  • NbOfRBVectors : number of RB mode vectors. Default: 6.

Definition at line 12 of file libmesh_solve_linear_system_input_parser.cpp.

13  {
14 
15  if (field_parser.search(1, "SysMatrix")) {
16  input_params.sys_matrix_file = field_parser.next(
17  input_params.sys_matrix_file);
18  } else {
19  homemade_error_msg("Missing the system matrix file!");
20  }
21 
22  if (field_parser.search(1, "SysRHSVector")) {
23  input_params.sys_rhs_vec_file = field_parser.next(
24  input_params.sys_rhs_vec_file);
25  } else {
26  homemade_error_msg("Missing the system RHS vector file!");
27  }
28 
29  if (field_parser.search(1, "OutputBase")) {
30  input_params.output_base = field_parser.next(
31  input_params.output_base);
32  } else {
33  homemade_error_msg("Missing the output filename base!");
34  }
35 
36  if (field_parser.search(1, "SysEps")) {
37  input_params.sys_eps = field_parser.next(
38  input_params.sys_eps);
39  } else {
40  input_params.sys_eps = 1e-5;
41  }
42 
43  if (field_parser.search(1, "SysIterDiv")) {
44  input_params.sys_iter_div = field_parser.next(
45  input_params.sys_iter_div);
46  } else {
47  input_params.sys_iter_div = 1000;
48  }
49 
50  if (field_parser.search(1, "RBVectorBase")) {
51  input_params.path_to_rb_vectors = field_parser.next(
52  input_params.path_to_rb_vectors);
53  input_params.bUseRBVectors = true;
54 
55  if (field_parser.search(1, "NbOfRBVectors")) {
56  input_params.nb_of_rb_vectors = field_parser.next(
57  input_params.nb_of_rb_vectors);
58  } else {
59  input_params.nb_of_rb_vectors = 6;
60  }
61 
62  } else {
63  input_params.bUseRBVectors = false;
64  }
65 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::get_input_params ( GetPot &  field_parser,
feti_setup_finish_params &  input_params 
)

Parser function for the coupled solver test programs.

Required parameters:

  • ClusterSchedulerType : scheduler type. Values: LOCAL, PBS or SLURM (code not implemented for the later yet, LOCAL runs the code without a scheduler).
  • ScratchFolderPath : path to the folder where the temporary files used by the coupled solver will be saved.
  • CouplingMatricesFolder : path to the folder containing the coupling matrices.

FETI / CG optional parameters:

  • CGPreconditionerType : CG preconditioner type. Values: "NONE", "Coupling_operator" or "Coupling_operator_jacobi".

Boolean flags:

  • UseRigidBodyModesB : use the rigid body modes for system B.

Rigid body mode parameters (only read if UseRigidBodyModesB is used):

  • RBVectorBase : filename base of the rigid body modes vectors.
  • NbOfRBVectors : number of RB mode vectors. Default: 6.

Definition at line 13 of file carl_feti_setup_finish_input_parser.cpp.

14  {
15 
16  if (field_parser.search(1, "ClusterSchedulerType")) {
17  std::string cluster_scheduler_type;
18  cluster_scheduler_type = field_parser.next(cluster_scheduler_type);
19  if(cluster_scheduler_type == "LOCAL")
20  {
21  std::cout << " !!! WARNING: " << std::endl;
22  std::cout << " Using the LOCAL job 'scheduler'. You will have to launch each script" << std::endl;
23  std::cout << " MANUALLY!!! Reason: MPI does not support recursive 'mpirun' calls" << std::endl;
24  input_params.scheduler = carl::ClusterSchedulerType::LOCAL;
25  }
26  else if(cluster_scheduler_type == "PBS")
27  input_params.scheduler = carl::ClusterSchedulerType::PBS;
28  else if(cluster_scheduler_type == "SLURM")
29  input_params.scheduler = carl::ClusterSchedulerType::SLURM;
30  else
31  homemade_error_msg("Invalid scheduler type!");
32  } else {
33  homemade_error_msg("Missing the scheduler type!");
34  }
35 
36  if (field_parser.search(1, "ScratchFolderPath")) {
37  input_params.scratch_folder_path = field_parser.next(
38  input_params.scratch_folder_path);
39  } else {
40  homemade_error_msg("Missing the external scratch folder path!");
41  }
42 
43  if (field_parser.search(1, "CouplingMatricesFolder")) {
44  input_params.coupling_folder_path = field_parser.next(
45  input_params.coupling_folder_path);
46  } else {
47  homemade_error_msg("Missing the coupling matrices path!");
48  }
49 
50  if (field_parser.search(1,"UseRigidBodyModesB"))
51  {
52  input_params.bUseRigidBodyModes = true;
53  if (field_parser.search(1, "RBVectorBase")) {
54  input_params.RB_vectors_base = field_parser.next(
55  input_params.RB_vectors_base);
56  } else {
57  homemade_error_msg("Missing the system B's rigid body mode vectors!");
58  }
59 
60  if (field_parser.search(1, "NbOfRBVectors")) {
61  input_params.nb_of_rb_vectors = field_parser.next(
62  input_params.nb_of_rb_vectors);
63  } else {
64  input_params.nb_of_rb_vectors = 6;
65  }
66  } else {
67  input_params.bUseRigidBodyModes = false;
68  }
69 
70  if ( field_parser.search(1, "CGPreconditionerType") )
71  {
72  std::string CG_precond_type_string = field_parser.next(CG_precond_type_string);
73  if(CG_precond_type_string == "NONE")
74  input_params.CG_precond_type = carl::BaseCGPrecondType::NO_PRECONDITIONER;
75  else if(CG_precond_type_string == "Coupling_operator")
76  input_params.CG_precond_type = carl::BaseCGPrecondType::COUPLING_OPERATOR;
77  else if(CG_precond_type_string == "Coupling_operator_jacobi")
78  input_params.CG_precond_type = carl::BaseCGPrecondType::COUPLING_JACOBI;
79  else
80  homemade_error_msg("Invalid preconditionner type!");
81  } else {
82  homemade_error_msg("Missing preconditionner type!");
83  }
84 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::get_input_params ( GetPot &  field_parser,
feti_setup_init_params &  input_params 
)

Parser function for the coupled solver test programs.

Required parameters:

  • Scheduler
    • ClusterSchedulerType : scheduler type. Values: LOCAL, PBS or SLURM (code not implemented for the later yet, LOCAL runs the code without a scheduler).
  • External solvers
    • ExtSolverA : command line for the external solver for system A.
    • ExtSolverB : command line for the external solver for system B.
    • ExtSolverAType : type of external solver used for system A. Values: LIBMESH_LINEAR.
    • ExtSolverBType : type of external solver used for system B. Values: LIBMESH_LINEAR.
    • ExtSolverAInput : path to a file containing the input parameters of the external solve for system A.
    • ExtSolverBInput : path to a file containing the input parameters of the external solve for system B.
  • Coupled solver
    • ScratchFolderPath : path to the folder where the temporary files used by the coupled solver will be saved.
    • CouplingMatricesFolder : path to the folder containing the coupling matrices.
    • OutputFolder : path to the coupled solution folder.

Boolean flags:

  • UseRigidBodyModesB : use the rigid body modes for system B.

Rigid body mode parameters (only read if UseRigidBodyModesB is used):

  • ExtForceSystemB : path to the vector containing the external forces for the system B.
  • RBVectorBase : filename base of the rigid body modes vectors.
  • NbOfRBVectors : number of RB mode vectors. Default: 6.

Optional parameters:

  • Scheduler
    • ScriptFile : path to the file used to generate the scripts, required if PBS or SLURM is used for the scheduler.
  • FETI / CG optional parameters:
    • CGPreconditionerType : CG preconditioner type. Values: "NONE", "Coupling_operator" or "Coupling_operator_jacobi". Default: "Coupling_operator".
    • CoupledConvAbs : CG absolute convergence on the residual. Default: 1e-20.
    • CoupledConvRel : CG relative convergence on the residual. Default: 1e-5.
    • CoupledCorrConvRel : CG relative convergence on the rigid body corrections. Default: 1e-6.
    • CoupledDiv : CG residual divergence parameter. Default: 100000.
    • CoupledIterMax : CG maximum number of iterations. Default: 1000.

Definition at line 13 of file carl_feti_setup_init_input_parser.cpp.

14  {
15 
16  if (field_parser.search(1, "ClusterSchedulerType")) {
17  std::string cluster_scheduler_type;
18  cluster_scheduler_type = field_parser.next(cluster_scheduler_type);
19  if(cluster_scheduler_type == "LOCAL")
20  {
21  std::cout << " !!! WARNING: " << std::endl;
22  std::cout << " The LOCAL 'scheduler' type is only intended for small and fast test cases" << std::endl;
23  std::cout << " on computers without a job scheduler (PBS, SLURM). You will have to launch" << std::endl;
24  std::cout << " each script MANUALLY!!! Reason: MPI does not support recursive 'mpirun' calls" << std::endl;
25  input_params.scheduler = carl::ClusterSchedulerType::LOCAL;
26  input_params.script_filename = "";
27  }
28  else if(cluster_scheduler_type == "PBS") {
29  input_params.scheduler = carl::ClusterSchedulerType::PBS;
30  if (field_parser.search(1, "ScriptFile")) {
31  input_params.script_filename = field_parser.next(input_params.script_filename);
32  } else {
33  homemade_error_msg("Missing the script file (needed for the PBS scheduler)!");
34  }
35  }
36  else if(cluster_scheduler_type == "SLURM") {
37  input_params.scheduler = carl::ClusterSchedulerType::SLURM;
38  if (field_parser.search(1, "ScriptFile")) {
39  input_params.script_filename = field_parser.next(input_params.script_filename);
40  } else {
41  homemade_error_msg("Missing the script file (needed for the SLURM scheduler)!");
42  }
43  }
44  else
45  homemade_error_msg("Invalid scheduler type!");
46  } else {
47  homemade_error_msg("Missing the scheduler type!");
48  }
49 
50  if (field_parser.search(1, "ExtSolverA")) {
51  input_params.ext_solver_BIG = field_parser.next(
52  input_params.ext_solver_BIG);
53  std::cout << input_params.ext_solver_BIG << std::endl;
54  } else {
55  homemade_error_msg("Missing the external solver A command line!");
56  }
57 
58  if (field_parser.search(1, "ExtSolverB")) {
59  input_params.ext_solver_micro = field_parser.next(
60  input_params.ext_solver_micro);
61  std::cout << input_params.ext_solver_micro << std::endl;
62  } else {
63  homemade_error_msg("Missing the external solver B command line!");
64  }
65 
66  std::string ext_solver_type;
67  if (field_parser.search(1, "ExtSolverAType")) {
68  ext_solver_type = field_parser.next(
69  ext_solver_type);
70  std::cout << ext_solver_type << std::endl;
71  if(ext_solver_type == "LIBMESH_LINEAR")
72  {
73  input_params.ext_solver_BIG_type = carl::ExtSolverType::LIBMESH_LINEAR;
74  } else {
75  homemade_error_msg("Invalid external solver A type!");
76  }
77  } else {
78  homemade_error_msg("Missing the external solver A type!");
79  }
80 
81  if (field_parser.search(1, "ExtSolverBType")) {
82  ext_solver_type = field_parser.next(
83  ext_solver_type);
84  std::cout << ext_solver_type << std::endl;
85  if(ext_solver_type == "LIBMESH_LINEAR")
86  {
87  input_params.ext_solver_BIG_type = carl::ExtSolverType::LIBMESH_LINEAR;
88  } else {
89  homemade_error_msg("Invalid external solver B type!");
90  }
91  } else {
92  homemade_error_msg("Missing the external solver B type!");
93  }
94 
95  if (field_parser.search(1, "ExtSolverAInput")) {
96  input_params.ext_solver_BIG_input = field_parser.next(
97  input_params.ext_solver_BIG_input);
98  std::cout << input_params.ext_solver_BIG_input << std::endl;
99  } else {
100  homemade_error_msg("Missing the input file for the external solver A!");
101  }
102 
103  if (field_parser.search(1, "ExtSolverBInput")) {
104  input_params.ext_solver_micro_input = field_parser.next(
105  input_params.ext_solver_micro_input);
106  std::cout << input_params.ext_solver_micro_input << std::endl;
107  } else {
108  homemade_error_msg("Missing the input file for the external solver B!");
109  }
110 
111  if (field_parser.search(1, "ScratchFolderPath")) {
112  input_params.scratch_folder_path = field_parser.next(
113  input_params.scratch_folder_path);
114  std::cout << input_params.scratch_folder_path << std::endl;
115  } else {
116  homemade_error_msg("Missing the external scratch folder path!");
117  }
118 
119  if (field_parser.search(1, "CouplingMatricesFolder")) {
120  input_params.coupling_folder_path = field_parser.next(
121  input_params.coupling_folder_path);
122  std::cout << input_params.coupling_folder_path << std::endl;
123  } else {
124  homemade_error_msg("Missing the coupling matrices path!");
125  }
126 
127  if (field_parser.search(1,"UseRigidBodyModesB"))
128  {
129  input_params.bUseRigidBodyModes = true;
130  if (field_parser.search(1, "ExtForceSystemB")) {
131  input_params.force_micro_path = field_parser.next(
132  input_params.force_micro_path);
133  std::cout << input_params.force_micro_path << std::endl;
134  } else {
135  homemade_error_msg("Missing the system B's external force file path!");
136  }
137 
138  if (field_parser.search(1, "RBVectorBase")) {
139  input_params.RB_vectors_base = field_parser.next(
140  input_params.RB_vectors_base);
141  std::cout << input_params.RB_vectors_base << std::endl;
142  } else {
143  homemade_error_msg("Missing the system B's rigid body mode vectors!");
144  }
145 
146  if (field_parser.search(1, "NbOfRBVectors")) {
147  input_params.nb_of_rb_vectors = field_parser.next(
148  input_params.nb_of_rb_vectors);
149  } else {
150  input_params.nb_of_rb_vectors = 6;
151  }
152  }
153  else
154  {
155  input_params.bUseRigidBodyModes = false;
156  }
157 
158  // Set CG coupling solver convergence
159  input_params.CG_coupled_conv_abs = 1e-20;
160  input_params.CG_coupled_conv_rel = 1e-5;
161  input_params.CG_coupled_div = 1e5;
162  input_params.CG_coupled_conv_max = 1e4;
163  input_params.CG_coupled_conv_corr =1e-5;
164 
165  if( field_parser.search(1,"CoupledConvAbs") )
166  {
167  input_params.CG_coupled_conv_abs = field_parser.next(input_params.CG_coupled_conv_abs);
168  }
169  if( field_parser.search(1,"CoupledConvRel") )
170  {
171  input_params.CG_coupled_conv_rel = field_parser.next(input_params.CG_coupled_conv_rel);
172  }
173  if( field_parser.search(1,"CoupledCorrConvRel") )
174  {
175  input_params.CG_coupled_conv_corr = field_parser.next(input_params.CG_coupled_conv_corr);
176  }
177  if( field_parser.search(1,"CoupledDiv") )
178  {
179  input_params.CG_coupled_div = field_parser.next(input_params.CG_coupled_div);
180  }
181  if( field_parser.search(1,"CoupledIterMax") )
182  {
183  input_params.CG_coupled_conv_max = field_parser.next(input_params.CG_coupled_conv_max);
184  }
185 
186  input_params.CG_precond_type = carl::BaseCGPrecondType::COUPLING_OPERATOR;
187  if ( field_parser.search(1, "CGPreconditionerType") )
188  {
189  std::string CG_precond_type_string = field_parser.next(CG_precond_type_string);
190  if(CG_precond_type_string == "NONE")
191  input_params.CG_precond_type = carl::BaseCGPrecondType::NO_PRECONDITIONER;
192  else if(CG_precond_type_string == "Coupling_operator")
193  input_params.CG_precond_type = carl::BaseCGPrecondType::COUPLING_OPERATOR;
194  else if(CG_precond_type_string == "Coupling_operator_jacobi")
195  input_params.CG_precond_type = carl::BaseCGPrecondType::COUPLING_JACOBI;
196  }
197 
198  if (field_parser.search(1,"OutputFolder")) {
199  input_params.output_folder = field_parser.next(
200  input_params.output_folder);
201  } else {
202  homemade_error_msg("Missing the output filename base!");
203  }
204 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::get_input_params ( GetPot &  field_parser,
feti_iterate_params &  input_params 
)

Parser function for the coupled solver test programs.

Required parameters:

  • ClusterSchedulerType : scheduler type. Values: LOCAL, PBS or SLURM (code not implemented for the later yet, LOCAL runs the code without a scheduler).
  • ScratchFolderPath : path to the folder where the temporary files used by the coupled solver will be saved.
  • CouplingMatricesFolder : path to the folder containing the coupling matrices.

Boolean flags:

  • UseRigidBodyModesB : use the rigid body modes for system B.

Rigid body mode parameters (only read if UseRigidBodyModesB is used):

  • ExtForceSystemB : path to the vector containing the external forces for the system B.
  • RBVectorBase : filename base of the rigid body modes vectors.

Optional parameters:

  • FETI / CG optional parameters:
    • CGPreconditionerType : CG preconditioner type. Values: "NONE", "Coupling_operator" or "Coupling_operator_jacobi". Default: "Coupling_operator".
    • CoupledConvAbs : CG absolute convergence on the residual. Default: 1e-20.
    • CoupledConvRel : CG relative convergence on the residual. Default: 1e-5.
    • CoupledCorrConvRel : CG relative convergence on the rigid body corrections. Default: 1e-6.
    • CoupledDiv : CG residual divergence parameter. Default: 100000.
    • CoupledIterMax : CG maximum number of iterations. Default: 1000.

Definition at line 13 of file carl_feti_iterate_input_parser.cpp.

14  {
15 
16  if (field_parser.search(1, "ClusterSchedulerType")) {
17  std::string cluster_scheduler_type;
18  cluster_scheduler_type = field_parser.next(cluster_scheduler_type);
19  if(cluster_scheduler_type == "LOCAL")
20  {
21  std::cout << " !!! WARNING: " << std::endl;
22  std::cout << " Using the LOCAL job 'scheduler'. You will have to launch each script" << std::endl;
23  std::cout << " MANUALLY!!! Reason: MPI does not support recursive 'mpirun' calls" << std::endl;
24  input_params.scheduler = carl::ClusterSchedulerType::LOCAL;
25  }
26  else if(cluster_scheduler_type == "PBS")
27  input_params.scheduler = carl::ClusterSchedulerType::PBS;
28  else if(cluster_scheduler_type == "SLURM")
29  input_params.scheduler = carl::ClusterSchedulerType::SLURM;
30  else
31  homemade_error_msg("Invalid scheduler type!");
32  } else {
33  homemade_error_msg("Missing the scheduler type!");
34  }
35 
36  if (field_parser.search(1, "ScratchFolderPath")) {
37  input_params.scratch_folder_path = field_parser.next(
38  input_params.scratch_folder_path);
39  } else {
40  homemade_error_msg("Missing the external scratch folder path!");
41  }
42 
43  if (field_parser.search(1, "CouplingMatricesFolder")) {
44  input_params.coupling_folder_path = field_parser.next(
45  input_params.coupling_folder_path);
46  } else {
47  homemade_error_msg("Missing the coupling matrices path!");
48  }
49 
50  if (field_parser.search(1,"UseRigidBodyModesB"))
51  {
52  input_params.bUseRigidBodyModes = true;
53  if (field_parser.search(1, "RBVectorBase")) {
54  input_params.RB_vectors_base = field_parser.next(
55  input_params.RB_vectors_base);
56  } else {
57  homemade_error_msg("Missing the system B's rigid body mode vectors!");
58  }
59 
60  if (field_parser.search(1, "NbOfRBVectors")) {
61  input_params.nb_of_rb_vectors = field_parser.next(
62  input_params.nb_of_rb_vectors);
63  } else {
64  input_params.nb_of_rb_vectors = 6;
65  }
66  }
67  else
68  {
69  input_params.bUseRigidBodyModes = false;
70  }
71 
72  // Set CG coupling solver convergence
73  input_params.CG_coupled_conv_abs = 1e-20;
74  input_params.CG_coupled_conv_rel = 1e-5;
75  input_params.CG_coupled_div = 1e5;
76  input_params.CG_coupled_conv_max = 1e4;
77  input_params.CG_coupled_conv_corr =1e-5;
78 
79  if( field_parser.search(1,"CoupledConvAbs") )
80  {
81  input_params.CG_coupled_conv_abs = field_parser.next(input_params.CG_coupled_conv_abs);
82  }
83  if( field_parser.search(1,"CoupledConvRel") )
84  {
85  input_params.CG_coupled_conv_rel = field_parser.next(input_params.CG_coupled_conv_rel);
86  }
87  if( field_parser.search(1,"CoupledCorrConvRel") )
88  {
89  input_params.CG_coupled_conv_corr = field_parser.next(input_params.CG_coupled_conv_corr);
90  }
91  if( field_parser.search(1,"CoupledDiv") )
92  {
93  input_params.CG_coupled_div = field_parser.next(input_params.CG_coupled_div);
94  }
95  if( field_parser.search(1,"CoupledIterMax") )
96  {
97  input_params.CG_coupled_conv_max = field_parser.next(input_params.CG_coupled_conv_max);
98  }
99 
100  if ( field_parser.search(1, "CGPreconditionerType") )
101  {
102  std::string CG_precond_type_string = field_parser.next(CG_precond_type_string);
103  if(CG_precond_type_string == "NONE")
104  input_params.CG_precond_type = carl::BaseCGPrecondType::NO_PRECONDITIONER;
105  else if(CG_precond_type_string == "Coupling_operator")
106  input_params.CG_precond_type = carl::BaseCGPrecondType::COUPLING_OPERATOR;
107  else if(CG_precond_type_string == "Coupling_operator_jacobi")
108  input_params.CG_precond_type = carl::BaseCGPrecondType::COUPLING_JACOBI;
109  else
110  homemade_error_msg("Invalid preconditionner type!");
111  } else {
112  homemade_error_msg("Missing preconditionner type!");
113  }
114 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::get_input_params ( GetPot &  field_parser,
feti_set_sol_params &  input_params 
)

Parser function for the coupled solver test programs.

Required parameters:

  • ScratchFolderPath : path to the folder where the temporary files used by the coupled solver will be saved.
  • OutputFolder : path to the coupled outup folder.

Boolean flags:

  • UseRigidBodyModesB : use the rigid body modes for system B.

Definition at line 14 of file carl_feti_set_sol_input_parser.cpp.

15  {
16 
17  if (field_parser.search(1, "ScratchFolderPath")) {
18  input_params.scratch_folder_path = field_parser.next(
19  input_params.scratch_folder_path);
20  } else {
21  homemade_error_msg("Missing the external scratch folder path!");
22  }
23 
24  if (field_parser.search(1,"UseRigidBodyModesB"))
25  {
26  input_params.bUseRigidBodyModes = true;
27  } else {
28  input_params.bUseRigidBodyModes = false;
29  }
30 
31  if (field_parser.search(1,"OutputFolder")) {
32  input_params.output_folder = field_parser.next(
33  input_params.output_folder);
34  } else {
35  homemade_error_msg("Missing the output filename base!");
36  }
37 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::get_input_params ( GetPot &  field_parser,
libmesh_apply_solution_input_params &  input_params 
)

Parser function for mesh deformation (apply solution) programs.

Required parameters:

  • InputVector or --input-vec: path to the input vector containing the displacements.
  • InputMesh or --input-mesh: path to the mesh that will be deformed.
  • PhysicalParameters or --physical-params: physical parameters file.
  • OutputMesh or --output-mesh: path to the output mesh.

Definition at line 12 of file libmesh_apply_solution_input_parser.cpp.

13  {
14 
15  if (field_parser.search(2, "InputVector", "--input-vec" )) {
16  input_params.input_vector = field_parser.next(
17  input_params.input_vector);
18  } else {
19  homemade_error_msg("Missing the displacement field vector!");
20  }
21 
22  if (field_parser.search(2, "InputMesh", "--input-mesh")) {
23  input_params.input_mesh = field_parser.next(
24  input_params.input_mesh);
25  } else {
26  homemade_error_msg("Missing the input mesh!");
27  }
28 
29  if (field_parser.search(2, "PhysicalParameters", "--physical-params")) {
30  input_params.physical_params_file = field_parser.next(
31  input_params.physical_params_file);
32  } else {
33  homemade_error_msg("Missing the physical params file!");
34  }
35 
36  if (field_parser.search(2, "OutputMesh", "--output-mesh")) {
37  input_params.output_mesh = field_parser.next(
38  input_params.output_mesh);
39  } else {
40  homemade_error_msg("Missing the output mesh path!");
41  }
42 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::get_intersection_input_params ( GetPot &  field_parser,
parallel_intersection_params &  input_params 
)

Parser function for the parallel intersection search program (source: CArl_build_intersections.cpp)

Required parameters:

  • MeshA, -mA or --meshA : path to the mesh A.
  • MeshB, -mB or --meshB : path to the mesh B.
  • MeshC, -mC or --meshC : path to the coupling mesh C.

Optional parameters:

  • OutputBase, -mO or --output : base of the output files (including folders). Default: test_inter.
  • MeshingMethod or --meshingMethodType : intersection meshing method. Values: CGAL or LIBMESH_TETGEN. Default: CGAL.

Boolean flags:

  • StitchInterMeshes : do not stich together the intersection meshes.
  • VerboseOutput or --verbose : print some extra information, such as the coupling mesh partitioning.

Definition at line 12 of file intersection_input_parser.cpp.

13  {
14 
15  // Set mesh files
16  if (field_parser.search(3, "--meshA", "-mA", "MeshA")) {
17  input_params.mesh_A = field_parser.next(
18  input_params.mesh_A);
19  } else {
20  homemade_error_msg("Missing mesh A!");
21  }
22 
23  if (field_parser.search(3, "--meshB", "-mB", "MeshB")) {
24  input_params.mesh_B = field_parser.next(
25  input_params.mesh_B);
26  } else {
27  homemade_error_msg("Missing mesh B!");
28  }
29 
30  if (field_parser.search(3, "--meshC", "-mC", "MeshC")) {
31  input_params.mesh_C = field_parser.next(
32  input_params.mesh_C);
33  } else {
34  homemade_error_msg("Missing the coupling mesh C!");
35  }
36 
37  if (field_parser.search(3, "--output", "-mO", "OutputBase")) {
38  input_params.output_base = field_parser.next(
39  input_params.output_base);
40  } else {
41  input_params.output_base = "inter_search";
42  }
43 
44  std::string meshing_method;
45  input_params.inter_meshing_method = carl::IntersectionMeshingMethod::CGAL;
46  if (field_parser.search(2, "--meshingMethodType", "MeshingMethod")) {
47  meshing_method = field_parser.next(
48  meshing_method);
49  if(meshing_method == "CGAL")
50  {
51  input_params.inter_meshing_method = carl::IntersectionMeshingMethod::CGAL;
52  }
53  else if(meshing_method == "TETGEN" )
54  {
55  input_params.inter_meshing_method = carl::IntersectionMeshingMethod::LIBMESH_TETGEN;
56  }
57  }
58 
59  if(field_parser.search(1,"StitchInterMeshes")) {
60  input_params.bStitchInterMeshes = true;
61  }
62  else
63  {
64  input_params.bStitchInterMeshes = false;
65  }
66 
67  if(field_parser.search(2,"VerboseOutput", "--verbose")) {
68  input_params.bVerbose = true;
69  }
70  else
71  {
72  input_params.bVerbose = false;
73  }
74 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::invert_index_unordered_map ( const std::unordered_map< int, int > &  input_map,
std::unordered_map< int, int > &  output_map 
)

Definition at line 72 of file common_functions.cpp.

75 {
76  int map_length = input_map.size();
77  output_map.reserve(map_length);
78 
79  std::unordered_map<int,int>::const_iterator mapIt = input_map.begin();
80  std::unordered_map<int,int>::const_iterator end_mapIt = input_map.end();
81 
82  for( ; mapIt != end_mapIt; ++mapIt)
83  {
84  output_map[mapIt->second] = mapIt->first;
85  }
86 }
template<typename T >
void carl::jump_lines ( T &  filestream,
unsigned int  numberOfLines = 1 
)

Definition at line 33 of file common_functions.h.

34 {
35  std::string dummy;
36  for(unsigned int iii = 0; iii < numberOfLines; ++iii)
37  std::getline(filestream,dummy);
38 };
void carl::lump_matrix ( libMesh::PetscMatrix< libMesh::Number > &  matrixInput,
libMesh::PetscMatrix< libMesh::Number > &  matrixOutput 
)

Definition at line 3 of file PETSC_matrix_operations.cpp.

5 {
6  if(matrixOutput.initialized())
7  {
8  matrixOutput.clear();
9  }
10 
11  int M = matrixInput.m();
12  int N = matrixInput.n();
13 
14  homemade_assert_msg(M == N, "Lumping: the matrix must be a square matrix");
15 
16  PetscInt local_M, local_N;
17 
18  MatGetLocalSize(matrixInput.mat(),&local_M,&local_N);
19 
20  // It will be a diagonal matrix, so no need of a heavy preallocation
21  matrixOutput.init(M,N,local_M,local_N,1,0);
22 
23  libMesh::PetscVector<libMesh::Number> UnityVec(matrixInput.comm(),M,local_M);
24  libMesh::PetscVector<libMesh::Number> DummyVector(matrixInput.comm(),M,local_M);
25 
26  VecSet(UnityVec.vec(),1);
27 
28  UnityVec.close();
29 
30  matrixInput.vector_mult(DummyVector,UnityVec);
31 
32  MatDiagonalSet(matrixOutput.mat(),DummyVector.vec(),INSERT_VALUES);
33 
34  matrixOutput.close();
35 }
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
void carl::lump_matrix_and_invert ( libMesh::PetscMatrix< libMesh::Number > &  matrixInput,
libMesh::PetscMatrix< libMesh::Number > &  matrixOutput 
)

Definition at line 37 of file PETSC_matrix_operations.cpp.

39 {
40  if(matrixOutput.initialized())
41  {
42  matrixOutput.clear();
43  }
44 
45  int M = matrixInput.m();
46  int N = matrixInput.n();
47 
48  homemade_assert_msg(M == N, "Lumping: the matrix must be a square matrix");
49 
50  PetscInt local_M, local_N;
51 
52  MatGetLocalSize(matrixInput.mat(),&local_M,&local_N);
53 
54  // It will be a diagonal matrix, so no need of a heavy preallocation
55  matrixOutput.init(M,N,local_M,local_N,1,0);
56 
57  libMesh::PetscVector<libMesh::Number> UnityVec(matrixInput.comm(),M,local_M);
58  libMesh::PetscVector<libMesh::Number> DummyVector(matrixInput.comm(),M,local_M);
59 
60  VecSet(UnityVec.vec(),1);
61 
62  UnityVec.close();
63 
64  matrixInput.vector_mult(DummyVector,UnityVec);
65 
66  DummyVector.reciprocal();
67 
68  MatDiagonalSet(matrixOutput.mat(),DummyVector.vec(),INSERT_VALUES);
69 
70  matrixOutput.close();
71 }
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
void carl::lump_matrix_and_invert ( libMesh::PetscMatrix< libMesh::Number > &  matrixInput,
libMesh::PetscVector< libMesh::Number > &  vecOutput 
)

Definition at line 73 of file PETSC_matrix_operations.cpp.

75 {
76  if(vecOutput.initialized())
77  {
78  vecOutput.clear();
79  }
80 
81  int M = matrixInput.m();
82  int N = matrixInput.n();
83 
84  homemade_assert_msg(M == N, "Lumping: the matrix must be a square matrix");
85 
86  PetscInt local_M, local_N;
87 
88  MatGetLocalSize(matrixInput.mat(),&local_M,&local_N);
89 
90  // It will be a diagonal matrix, so no need of a heavy preallocation
91  vecOutput.init(M,local_M,false);
92 
93  libMesh::PetscVector<libMesh::Number> UnityVec(matrixInput.comm(),M,local_M);
94 
95  VecSet(UnityVec.vec(),1);
96 
97  UnityVec.close();
98 
99  matrixInput.vector_mult(vecOutput,UnityVec);
100 
101  vecOutput.reciprocal();
102 }
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
template<typename T >
void carl::MPI_reduce_vector ( std::vector< T > &  r,
int  root,
const libMesh::Parallel::Communicator &  Comm 
)
inline

Definition at line 22 of file mpi_carl_tools.h.

23 {
24  if (Comm.size() > 1 && !r.empty())
25  {
26  libmesh_assert(Comm.verify(r.size()));
27 
28  std::vector<T> temp(r);
29  MPI_Reduce(&temp[0], &r[0], libMesh::cast_int<int>(r.size()),
30  libMesh::Parallel::StandardType<T>(&temp[0]), MPI_SUM, root,
31  Comm.get());
32  }
33  }
void carl::PETSC_invert_dense_matrix ( Mat &  matrix_in,
Mat &  matrix_out 
)

Definition at line 490 of file PETSC_matrix_operations.cpp.

491 {
492  // WARNING: this operation is extremely expensive for large systems !!!
493  // Only use this for small matrices (example: the 6x6 or 3x3
494  // or 3x3 matrices from the null space projectors,
495 
496  Mat Id_mat;
497  MatFactorInfo factor_info;
498  IS rperm, cperm;
499 
500  // Duplicate matrix data structures
501  MatDuplicate(matrix_in,MAT_DO_NOT_COPY_VALUES,&matrix_out);
502  MatDuplicate(matrix_in,MAT_DO_NOT_COPY_VALUES,&Id_mat);
503 
504  // Create identity
505  MatZeroEntries(Id_mat);
506  MatShift(Id_mat,1);
507 
508  // Factor input matrix
509  MatGetOrdering(matrix_in,MATORDERINGNATURAL, &rperm, &cperm);
510  MatFactorInfoInitialize(&factor_info);
511  MatLUFactor(matrix_in,rperm,cperm,&factor_info);
512 
513  // Calculate inverse
514  MatMatSolve(matrix_in,Id_mat,matrix_out);
515 
516  // Reset input's factoring
517  MatSetUnfactored(matrix_in);
518 
519  // Cleanup
520  MatDestroy(&Id_mat);
521 }
void carl::PETSC_MatMultScale_Bcast ( Mat  mat_seq,
Vec  vec_seq_in,
Vec  vec_seq_out,
double  a_const 
)

Definition at line 523 of file PETSC_matrix_operations.cpp.

524 {
525  // Calculate the product vec_seq_out = a_const * (mat_seq * vec_seq_in)
526  // on the first processor, and then sync the result
527 
528  // First, check if the matrices and vectors are sequential
529  MatType mat_type_query;
530  VecType vec_in_type_query;
531  VecType vec_out_type_query;
532 
533  MatGetType(mat_seq,&mat_type_query);
534  VecGetType(vec_seq_in,&vec_in_type_query);
535  VecGetType(vec_seq_out,&vec_out_type_query);
536 
537  homemade_assert_msg(std::strcmp(mat_type_query,MATSEQDENSE) == 0 ,"Matrix is not dense and sequential");
538  homemade_assert_msg(std::strcmp(vec_in_type_query,VECSEQ) == 0 ,"Input vector is not sequential");
539  homemade_assert_msg(std::strcmp(vec_out_type_query,VECSEQ) == 0 ,"Output vector is not sequential");
540 
541  // Get the vector dimension
542  PetscInt vec_dim = 0;
543  VecGetSize(vec_seq_in,&vec_dim);
544 
545  // Get the ranks
546  int rank;
547  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
548 
549  // Set the dummy pointer and vector
550  PetscScalar * dummy_seq_array;
551 
552  // Do the product in the first processor
553  if(rank == 0)
554  {
555  MatMult(mat_seq,vec_seq_in,vec_seq_out);
556  VecScale(vec_seq_out,a_const);
557  }
558 
559  // Sync
560  VecGetArray(vec_seq_out,&dummy_seq_array);
561  MPI_Bcast(dummy_seq_array, vec_dim, MPIU_SCALAR, 0, PETSC_COMM_WORLD);
562  MPI_Barrier(PETSC_COMM_WORLD);
563  VecRestoreArray(vec_seq_out,&dummy_seq_array);
564 }
homemade_assert_msg
#define homemade_assert_msg(asserted, msg)
Definition: common_header.h:63
void carl::print_input_params ( const std::string &  output_filename,
libmesh_solve_linear_system_input_params &  input_params 
)

Function used to generate a solver input file from "input_params".

Definition at line 67 of file libmesh_solve_linear_system_input_parser.cpp.

68  {
69 
70  std::ofstream output_file(output_filename);
71 
72  output_file << "SysMatrix " << input_params.sys_matrix_file << std::endl;
73  output_file << "SysRHSVector " << input_params.sys_rhs_vec_file << std::endl;
74  output_file << "OutputBase " << input_params.output_base << std::endl;
75 
76  output_file << "SysEps " << input_params.sys_eps << std::endl;
77  output_file << "SysIterDiv " << input_params.sys_iter_div << std::endl;
78 
79  if(input_params.bUseRBVectors)
80  {
81  output_file << "RBVectorBase " << input_params.path_to_rb_vectors << std::endl;
82  output_file << "NbOfRBVectors " << input_params.nb_of_rb_vectors << std::endl;
83  }
84 
85  output_file.close();
86 };
void carl::print_matrix ( libMesh::PetscMatrix< libMesh::Number > &  CouplingTestMatrix)

Definition at line 104 of file PETSC_matrix_operations.cpp.

105 {
106  libMesh::Real accumulator = 0;
107  const libMesh::Parallel::Communicator& MatrixComm = CouplingTestMatrix.comm();
108 
109  std::cout << "| M_i,j : " << CouplingTestMatrix.m() << " x " << CouplingTestMatrix.n() << std::endl;
110 
111  int nodes = MatrixComm.size();
112 
113  PetscInt local_M, local_N;
114  MatGetLocalSize(CouplingTestMatrix.mat(),&local_M,&local_N);
115 
116  bool bPrintOnTerminal = CouplingTestMatrix.m() < 101 && CouplingTestMatrix.n() < 101 && nodes > 1;
117  if(bPrintOnTerminal)
118  {
119  for(unsigned int iii = 0; iii < CouplingTestMatrix.m(); ++iii)
120  {
121  for(unsigned int jjj = 0; jjj < CouplingTestMatrix.n(); ++jjj)
122  {
123  std::cout << " " << CouplingTestMatrix(iii,jjj);
124  }
125  std::cout << std::endl;
126  }
127  }
128 
129  libMesh::PetscVector<libMesh::Number> dummy_vec(MatrixComm,CouplingTestMatrix.n(),local_N);
130  MatGetRowSum(CouplingTestMatrix.mat(),dummy_vec.vec());
131 
132  VecSum(dummy_vec.vec(),&accumulator);
133  std::cout << "|" << std::endl;
134  std::cout << "| Sum( M_i,j ) = " << accumulator << std::endl << std::endl;
135 }
void carl::print_matrix_col_line_sum ( libMesh::PetscMatrix< libMesh::Number > &  CouplingTestMatrix,
const std::string  name_base 
)

Definition at line 137 of file PETSC_matrix_operations.cpp.

138 {
139  libMesh::PetscVector<libMesh::Number> col_sum(CouplingTestMatrix.comm(),CouplingTestMatrix.m());
140  libMesh::PetscVector<libMesh::Number> row_sum(CouplingTestMatrix.comm(),CouplingTestMatrix.n());
141 
142  for(unsigned int iii = 0; iii < CouplingTestMatrix.m(); ++iii)
143  {
144  for(unsigned int jjj = 0; jjj < CouplingTestMatrix.n(); ++jjj)
145  {
146  col_sum.add(iii,CouplingTestMatrix(iii,jjj));
147  row_sum.add(jjj,CouplingTestMatrix(iii,jjj));
148  }
149  }
150 
151 // Print MatLab debugging output? Variable defined at "carl_headers.h"
152 #ifdef PRINT_MATLAB_DEBUG
153  col_sum.print_matlab(name_base + "_col.m");
154  row_sum.print_matlab(name_base + "_row.m");
155 #endif
156 }
void carl::print_matrix_dim ( libMesh::PetscMatrix< libMesh::Number > &  CouplingTestMatrix,
bool  bDetailed = false 
)

Definition at line 165 of file PETSC_matrix_operations.cpp.

166 {
167  std::cout << "| M_i,j : " << CouplingTestMatrix.m() << " x " << CouplingTestMatrix.n() << std::endl;
168 // if(bDetailed)
169 // {
170  MatInfo temp_info;
171  MatGetInfo(CouplingTestMatrix.mat(),MAT_LOCAL,&temp_info);
172  std::cout << "| LOCAL : memory = " << temp_info.memory << std::endl;
173  std::cout << "| non-zeros used = " << (100.*temp_info.nz_used)/temp_info.nz_allocated << " % " << std::endl;
174 
175 // int non_zeros = temp_info.nz_used;
176 // std::vector<int> all_temp_info;
177 // CouplingTestMatrix.comm().gather(0,non_zeros,all_temp_info);
178 
179 // if(CouplingTestMatrix.comm().rank() == 0)
180 // {
181 // std::cout << " | ";
182 // for(unsigned int iii = 0; iii < CouplingTestMatrix.comm().size(); ++iii)
183 // {
184 // std::cout << all_temp_info[iii] << " ";
185 // }
186 // std::cout << std::endl << std::endl;
187 // }
188 
189  MatGetInfo(CouplingTestMatrix.mat(),MAT_GLOBAL_SUM,&temp_info);
190  std::cout << "| GLOBAL : memory = " << temp_info.memory << std::endl;
191  std::cout << "| non-zeros used = " << (100.*temp_info.nz_used)/temp_info.nz_allocated << " % " << std::endl;
192  std::cout << "| total nb. of non-zeros used = " << temp_info.nz_used << std::endl;
193  std::cout << "| total nb. of non-zeros allocated = " << temp_info.nz_allocated << std::endl;
194  std::cout << "| sparsity = " << (100.*temp_info.nz_used)/(CouplingTestMatrix.m() * CouplingTestMatrix.n()) << " % " << std::endl;
195 
196  MatGetInfo(CouplingTestMatrix.mat(),MAT_GLOBAL_MAX,&temp_info);
197  std::cout << "| MAX : memory = " << temp_info.memory << std::endl;
198  std::cout << "| non-zeros used = " << (100.*temp_info.nz_used)/temp_info.nz_allocated << " % " << std::endl << std::endl;
199  //}
200 }
void carl::print_matrix_info ( libMesh::PetscMatrix< libMesh::Number > &  InputMatrix,
std::ostream &  os = libMesh::out 
)

Definition at line 202 of file PETSC_matrix_operations.cpp.

203 {
204  const libMesh::Parallel::Communicator& WorldComm = InputMatrix.comm();
205 
206  unsigned int rank = WorldComm.rank();
207  unsigned int nodes = WorldComm.size();
208 
209  MatInfo local_info;
210  MatInfo global_info;
211  MatGetInfo(InputMatrix.mat(),MAT_LOCAL,&local_info);
212  MatGetInfo(InputMatrix.mat(),MAT_LOCAL,&global_info);
213 
214  // Set up local variables
215  int l_nz_used = local_info.nz_used;
216  int l_nz_allocated = local_info.nz_allocated;
217  int l_memory = local_info.memory;
218 
219  int l_n = -1;
220  int l_m = -1;
221 
222  double accumulator_n = 0;
223  double accumulator_m = 0;
224 
225  MatGetLocalSize(InputMatrix.mat(),&l_n,&l_m);
226 
227  std::vector<int> full_nz_used;
228  std::vector<int> full_nz_allocated;
229  std::vector<int> full_memory;
230 
231  std::vector<int> full_n;
232  std::vector<int> full_m;
233 
234  if(rank == 0)
235  {
236  full_nz_used.resize(nodes,0);
237  full_nz_allocated.resize(nodes,0);
238  full_memory.resize(nodes,0);
239  full_n.resize(nodes,0);
240  full_m.resize(nodes,0);
241  }
242 
243  WorldComm.gather(0,l_nz_used,full_nz_used);
244  WorldComm.gather(0,l_nz_allocated,full_nz_allocated);
245  WorldComm.gather(0,l_memory,full_memory);
246  WorldComm.gather(0,l_n,full_n);
247  WorldComm.gather(0,l_m,full_m);
248 
249  if(rank == 0)
250  {
251  os << "# rank \t local_n \t local_m \t start_n \t start_m \t nz_alloc \t nz_used \t memory " << std::endl;
252  for(unsigned int iii = 0; iii < nodes; ++iii)
253  {
254  os << iii << " \t "
255  << full_n[iii] << " \t "
256  << full_m[iii] << " \t "
257  << accumulator_n << " \t "
258  << accumulator_m << " \t "
259  << full_nz_allocated[iii] << " \t "
260  << full_nz_used[iii] << " \t "
261  << full_memory[iii] << std::endl;
262  accumulator_n += full_n[iii];
263  accumulator_m += full_m[iii];
264  }
265  }
266 
267  WorldComm.barrier();
268 }
void carl::print_matrix_matlab ( libMesh::PetscMatrix< libMesh::Number > &  CouplingTestMatrix,
const std::string  name_base 
)

Definition at line 158 of file PETSC_matrix_operations.cpp.

159 {
160  std::cout << "| M_i,j : " << CouplingTestMatrix.m() << " x " << CouplingTestMatrix.n() << std::endl;
161 
162  CouplingTestMatrix.print_matlab(name_base);
163 }
void carl::print_PETSC_vector ( libMesh::PetscVector< libMesh::Number > &  input_vec,
const std::string &  filename 
)
void carl::print_stats_to_file ( std::vector< double > &  vec_data,
const std::string  filename 
)

Definition at line 126 of file common_functions.cpp.

127 {
128  std::ofstream output_stream(filename,std::ofstream::app);
129  libMesh::StatisticsVector<double> statistics_vec(vec_data.size(),0);
130  for(unsigned int iii = 0; iii < vec_data.size(); ++iii)
131  {
132  statistics_vec[iii] = vec_data[iii];
133  }
134 
135  output_stream << statistics_vec.minimum() << " "
136  << statistics_vec.maximum() << " "
137  << statistics_vec.mean() << " "
138  << statistics_vec.median() << " "
139  << statistics_vec.stddev() << std::endl;
140 
141  output_stream.close();
142 };
void carl::read_local_intersection_tables ( const libMesh::Parallel::Communicator &  WorldComm,
const std::string &  intersection_local_table_Filename,
std::unordered_map< int, std::pair< int, int > > &  local_intersection_pairs_map,
std::unordered_map< int, int > &  local_intersection_meshI_to_inter_map 
)

Definition at line 546 of file mesh_tables.cpp.

552 {
553  int nbOfIntersections = -1;
554  int nbOfInterElems = -1;
555 
556  // Declare a few auxiliary variables
557  int temp_interID = -1;
558  int temp_idxA = -1;
559  int temp_idxB = -1;
560  int temp_idxI = -1;
561  int temp_nbOfInter = -1;
562  int temp_dummy = -1;
563 
564  int interIdx = 0;
565 
566  // Do the file reading
567  std::ifstream intersection_full_file(intersection_local_table_Filename);
568 
569  intersection_full_file >> nbOfIntersections >> nbOfInterElems >> temp_dummy;
570 
571  for(int iii = 0; iii < nbOfIntersections; ++iii)
572  {
573  // For each line, read ...
574 
575  // ... the inter ID, A and B's elements, and the number of I's elements ...
576  intersection_full_file >> temp_interID
577  >> temp_idxA >> temp_idxB
578  >> temp_nbOfInter;
579 
580  local_intersection_pairs_map[temp_interID] = std::pair<int,int>(temp_idxA,temp_idxB);
581 
582  // ... and all of I's elements
583  for(int jjj = 0; jjj < temp_nbOfInter; ++jjj)
584  {
585  intersection_full_file >> temp_idxI;
586  local_intersection_meshI_to_inter_map[temp_idxI] = temp_interID;
587  ++interIdx;
588  }
589  }
590  intersection_full_file.close();
591 
592  WorldComm.barrier();
593 };
void carl::read_PETSC_matrix ( Mat  input_mat,
const std::string &  filename,
MPI_Comm  comm = PETSC_COMM_WORLD 
)

Definition at line 439 of file PETSC_matrix_operations.cpp.

441 {
442 
443  PetscViewer viewer;
444  PetscViewerBinaryOpen(comm,filename.c_str(),FILE_MODE_READ,&viewer);
445  MatLoad(input_mat,viewer);
446 
447  PetscViewerDestroy(&viewer);
448 }
void carl::read_PETSC_matrix ( libMesh::PetscMatrix< libMesh::Number > &  input_mat,
const std::string &  filename 
)

Definition at line 450 of file PETSC_matrix_operations.cpp.

452 {
453  read_PETSC_matrix(input_mat.mat(),filename,input_mat.comm().get());
454 }
carl::read_PETSC_matrix
void read_PETSC_matrix(Mat input_mat, const std::string &filename, MPI_Comm comm=PETSC_COMM_WORLD)
Definition: PETSC_matrix_operations.cpp:439
void carl::read_PETSC_vector ( libMesh::PetscVector< libMesh::Number > &  input_vec,
const std::string &  filename 
)

Definition at line 377 of file PETSC_matrix_operations.cpp.

379 {
380  read_PETSC_vector(input_vec.vec(),filename,input_vec.comm().get());
381 };
carl::read_PETSC_vector
void read_PETSC_vector(libMesh::PetscVector< libMesh::Number > &input_vec, const std::string &filename)
Definition: PETSC_matrix_operations.cpp:377
void carl::read_PETSC_vector ( Vec  input_vec,
const std::string &  filename,
MPI_Comm  comm = PETSC_COMM_WORLD 
)

Definition at line 367 of file PETSC_matrix_operations.cpp.

369 {
370  PetscViewer viewer;
371  PetscViewerBinaryOpen(comm,filename.c_str(),FILE_MODE_READ,&viewer);
372  VecLoad(input_vec,viewer);
373 
374  PetscViewerDestroy(&viewer);
375 }
template<typename Sys >
void carl::reduced_system_init ( Sys &  system_input)

Definition at line 74 of file common_functions.h.

75 {
76  libMesh::DofMap& system_dof_map = system_input.get_dof_map();
77  libMesh::MeshBase& system_mesh = system_input.get_mesh();
78 
79  unsigned int nb_of_variable_groups = system_input.n_variable_groups();
80  for (unsigned int vg=0; vg<nb_of_variable_groups; vg++)
81  {
82  system_dof_map.add_variable_group(system_input.variable_group(vg));
83  }
84 
85  system_dof_map.distribute_dofs(system_mesh);
86  system_input.reinit_constraints();
87  system_dof_map.prepare_send_list();
88  system_dof_map.compute_sparsity(system_mesh);
89 };
void carl::repartition_system_meshes ( const libMesh::Parallel::Communicator &  WorldComm,
libMesh::Mesh &  mesh_input,
libMesh::Mesh &  mesh_intersect,
std::unordered_map< int, std::pair< int, int > > &  local_intersection_pairs_map,
bool  bUseSecond = true 
)

Definition at line 755 of file mesh_tables.cpp.

761 {
762  std::vector<double> mesh_input_weights;
763  mesh_input_weights.resize(mesh_input.n_elem(),0);
764 
765  // Nb. of intersections associated to each elem of mesh_intersect
766  std::vector<int> inter_per_elem(mesh_intersect.n_elem());
767 
768  std::unordered_map<int,std::pair<int,int> >::const_iterator inter_pairsIt =
769  local_intersection_pairs_map.begin();
770  const std::unordered_map<int,std::pair<int,int> >::const_iterator end_inter_pairsIt =
771  local_intersection_pairs_map.end();
772  for( ; inter_pairsIt != end_inter_pairsIt; ++inter_pairsIt)
773  {
774  const std::pair<int,int>& dummy_pair = inter_pairsIt->second;
775 
776  if(bUseSecond)
777  {
778  ++inter_per_elem[dummy_pair.first];
779  }
780  else
781  {
782  ++inter_per_elem[dummy_pair.second];
783  }
784  }
785 
786  inter_pairsIt = local_intersection_pairs_map.begin();
787  for( ; inter_pairsIt != end_inter_pairsIt; ++inter_pairsIt)
788  {
789  const std::pair<int,int>& dummy_pair = inter_pairsIt->second;
790 
791  if(bUseSecond)
792  {
793  mesh_input_weights[dummy_pair.second] += inter_per_elem[dummy_pair.first];
794  }
795  else
796  {
797  mesh_input_weights[dummy_pair.first] += inter_per_elem[dummy_pair.second];
798  }
799  }
800 
801  WorldComm.sum(mesh_input_weights);
802 
803  // Basic partitioning
804  libMesh::MeshBase::const_element_iterator mesh_inputIt =
805  mesh_input.elements_begin();
806  const libMesh::MeshBase::const_element_iterator end_mesh_inputIt =
807  mesh_input.elements_end();
808  for( ; mesh_inputIt != end_mesh_inputIt; ++mesh_inputIt )
809  {
810  const libMesh::Elem* elem_input = *mesh_inputIt;
811  int idxinput = elem_input->id();
812 
813  mesh_input_weights[idxinput] += elem_input->n_nodes();
814  }
815 
816  libMesh::ErrorVector dummy_input_error(mesh_input.n_elem());
817  for(unsigned int iii = 0; iii < mesh_input.n_elem(); ++iii)
818  dummy_input_error[iii] = mesh_input_weights[iii];
819 
820  libMesh::Partitioner * dummy_partitioner_input = mesh_input.partitioner().get();
821  dummy_partitioner_input->attach_weights(&dummy_input_error);
822  mesh_input.partition();
823 };
void carl::set_equivalence_tables ( const libMesh::Parallel::Communicator &  WorldComm,
const std::string &  equivalence_table_A_Filename,
const std::string &  equivalence_table_B_Filename,
std::unordered_map< int, int > &  equivalence_table_A_to_R_A,
std::unordered_map< int, int > &  equivalence_table_B_to_R_B,
std::unordered_map< int, int > &  equivalence_table_R_A_to_A,
std::unordered_map< int, int > &  equivalence_table_R_B_to_B 
)

Definition at line 306 of file mesh_tables.cpp.

315 {
316  int rank = WorldComm.rank();
317 
318  // While the equivalence tables are saved as unordered maps, it's easier to
319  // save them as vectors at first, broadcast them, and them reconvert to maps
320  std::vector<int> dummy_equivalence_table_A;
321  std::vector<int> dummy_equivalence_table_B;
322 
323  int nbOfRestricted_A_Elems = -1;
324  int nbOfRestricted_B_Elems = -1;
325 
326  // Read files with proc 0
327  if(rank == 0)
328  {
329  int temp_RX = -1;
330  int temp_X = -1;
331 
332  std::ifstream equivalence_A_file(equivalence_table_A_Filename);
333 
334  equivalence_A_file >> nbOfRestricted_A_Elems;
335  dummy_equivalence_table_A.resize(2*nbOfRestricted_A_Elems);
336 
337  for(int iii = 0; iii < nbOfRestricted_A_Elems; ++iii)
338  {
339  equivalence_A_file >> temp_RX
340  >> temp_X;
341 
342  dummy_equivalence_table_A[2*iii] = temp_RX;
343  dummy_equivalence_table_A[2*iii + 1] = temp_X;
344  }
345  equivalence_A_file.close();
346 
347  std::ifstream equivalence_B_file(equivalence_table_B_Filename);
348 
349  equivalence_B_file >> nbOfRestricted_B_Elems;
350  dummy_equivalence_table_B.resize(2*nbOfRestricted_B_Elems);
351 
352  for(int iii = 0; iii < nbOfRestricted_B_Elems; ++iii)
353  {
354  equivalence_B_file >> temp_RX
355  >> temp_X;
356 
357  dummy_equivalence_table_B[2*iii] = temp_RX;
358  dummy_equivalence_table_B[2*iii + 1] = temp_X;
359  }
360  equivalence_B_file.close();
361  }
362 
363  // Broadcast the sizes and resize on other procs
364  WorldComm.broadcast(nbOfRestricted_A_Elems);
365  WorldComm.broadcast(nbOfRestricted_B_Elems);
366 
367  if(rank != 0)
368  {
369  dummy_equivalence_table_A.resize(2*nbOfRestricted_A_Elems);
370  dummy_equivalence_table_B.resize(2*nbOfRestricted_B_Elems);
371  }
372 
373  WorldComm.barrier();
374 
375  WorldComm.broadcast(dummy_equivalence_table_A);
376  WorldComm.broadcast(dummy_equivalence_table_B);
377 
378  // Convert back to unoredered maps
379  equivalence_table_A_to_R_A.reserve(nbOfRestricted_A_Elems);
380  equivalence_table_B_to_R_B.reserve(nbOfRestricted_B_Elems);
381 
382  equivalence_table_R_A_to_A.reserve(nbOfRestricted_A_Elems);
383  equivalence_table_R_B_to_B.reserve(nbOfRestricted_B_Elems);
384 
385  for(int iii = 0; iii < nbOfRestricted_A_Elems; ++iii)
386  {
387  equivalence_table_R_A_to_A[dummy_equivalence_table_A[2*iii]] =
388  dummy_equivalence_table_A[2*iii + 1];
389  equivalence_table_A_to_R_A[dummy_equivalence_table_A[2*iii + 1]] =
390  dummy_equivalence_table_A[2*iii];
391  }
392 
393  for(int iii = 0; iii < nbOfRestricted_B_Elems; ++iii)
394  {
395  equivalence_table_R_B_to_B[dummy_equivalence_table_B[2*iii]] =
396  dummy_equivalence_table_B[2*iii + 1];
397  equivalence_table_B_to_R_B[dummy_equivalence_table_B[2*iii + 1]] =
398  dummy_equivalence_table_B[2*iii];
399  }
400 }
void carl::set_full_intersection_tables ( const libMesh::Parallel::Communicator &  WorldComm,
const std::string &  intersection_full_table_Filename,
std::unordered_map< int, std::pair< int, int > > &  full_intersection_pairs_map,
std::unordered_map< int, int > &  full_intersection_meshI_to_inter_map 
)

Definition at line 437 of file mesh_tables.cpp.

443 {
444  int rank = WorldComm.rank();
445 
446  // While the equivalence tables are saved as unordered maps, it's easier to
447  // save them as vectors at first, broadcast them, and them reconvert to maps
448  std::vector<int> dummy_intersections_IDs;
449  std::vector<int> dummy_intersection_pairs_table;
450  std::vector<int> dummy_all_intersections_table;
451  std::vector<int> dummy_intersections_sizes_table;
452 
453  int nbOfIntersections = -1;
454  int nbOfInterElems = -1;
455 
456  // Declare a few auxiliary variables
457  int temp_interID = -1;
458  int temp_idxA = -1;
459  int temp_idxB = -1;
460  int temp_idxI = -1;
461  int temp_nbOfInter = -1;
462 
463  int interIdx = 0;
464 
465  // Do the file reading job on proc 0
466  if(rank == 0)
467  {
468  std::ifstream intersection_full_file(intersection_full_table_Filename);
469 
470  intersection_full_file >> nbOfIntersections >> nbOfInterElems;
471 
472  dummy_intersections_IDs.resize(nbOfIntersections);
473  dummy_intersection_pairs_table.resize(2*nbOfIntersections);
474  dummy_all_intersections_table.resize(nbOfInterElems);
475  dummy_intersections_sizes_table.resize(nbOfIntersections);
476 
477  for(int iii = 0; iii < nbOfIntersections; ++iii)
478  {
479  // For each line, read ...
480 
481  // ... the inter ID, A and B's elements, and the number of I's elements ...
482  intersection_full_file >> temp_interID
483  >> temp_idxA >> temp_idxB
484  >> temp_nbOfInter;
485 
486  dummy_intersections_IDs[iii] = temp_interID;
487  dummy_intersection_pairs_table[2*iii] = temp_idxA;
488  dummy_intersection_pairs_table[2*iii + 1] = temp_idxB;
489  dummy_intersections_sizes_table[iii] = temp_nbOfInter;
490 
491  // ... and all of I's elements
492  for(int jjj = 0; jjj < temp_nbOfInter; ++jjj)
493  {
494  intersection_full_file >> temp_idxI;
495  dummy_all_intersections_table[interIdx] = temp_idxI;
496 
497  ++interIdx;
498  }
499  }
500  intersection_full_file.close();
501  }
502 
503  // Broadcast the sizes and resize on other procs
504  WorldComm.broadcast(nbOfIntersections);
505  WorldComm.broadcast(nbOfInterElems);
506 
507  if(rank != 0)
508  {
509  dummy_intersections_IDs.resize(nbOfIntersections);
510  dummy_intersection_pairs_table.resize(2*nbOfIntersections);
511  dummy_all_intersections_table.resize(nbOfInterElems);
512  dummy_intersections_sizes_table.resize(nbOfIntersections);
513  }
514 
515  WorldComm.barrier();
516 
517  WorldComm.broadcast(dummy_intersections_IDs);
518  WorldComm.broadcast(dummy_intersection_pairs_table);
519  WorldComm.broadcast(dummy_all_intersections_table);
520  WorldComm.broadcast(dummy_intersections_sizes_table);
521 
522  full_intersection_pairs_map.reserve(nbOfIntersections);
523  full_intersection_meshI_to_inter_map.reserve(nbOfInterElems);
524 
525  interIdx = 0;
526 
527  // Convert the vectors to maps
528  for(int iii = 0; iii < nbOfIntersections; ++iii)
529  {
530  temp_interID = dummy_intersections_IDs[iii];
531  temp_idxA = dummy_intersection_pairs_table[2*iii];
532  temp_idxB = dummy_intersection_pairs_table[2*iii + 1];
533 
534  full_intersection_pairs_map[temp_interID] = std::pair<int,int>(temp_idxA,temp_idxB);
535 
536  temp_nbOfInter = dummy_intersections_sizes_table[iii];
537  for(int jjj = 0; jjj < temp_nbOfInter; ++jjj)
538  {
539  temp_idxI = dummy_all_intersections_table[interIdx];
540  full_intersection_meshI_to_inter_map[temp_idxI] = temp_interID;
541  ++interIdx;
542  }
543  }
544 };
void carl::set_global_mediator_system_intersection_lists ( const libMesh::Parallel::Communicator &  WorldComm,
const std::string &  intersection_global_table_Filename,
const std::unordered_map< int, int > &  equivalence_table_system_to_mediator,
const std::unordered_map< int, int > &  equivalence_table_mediator_to_system,
std::unordered_multimap< int, int > &  inter_mediator_A,
std::unordered_multimap< int, int > &  inter_mediator_B 
)

Definition at line 683 of file mesh_tables.cpp.

691 {
692  // First, do the reading work on proc. 0
693  int rank = WorldComm.rank();
694  int nodes = WorldComm.size();
695 
696  int temp_interID;
697  std::vector<int> temp_idxA;
698  std::vector<int> temp_idxB;
699 
700  int nbOfIntersections = -1;
701  int nbOfInterElems = -1;
702 
703  std::string dummy_data;
704 
705  if(rank == 0)
706  {
707  std::ifstream intersection_full_file(intersection_global_table_Filename);
708  intersection_full_file >> nbOfIntersections >> nbOfInterElems;
709 
710  temp_idxA.resize(nbOfIntersections);
711  temp_idxB.resize(nbOfIntersections);
712 
713  for(int iii = 0; iii < nbOfIntersections; ++iii)
714  {
715  // intersection_full_file >> temp_interID
716  // >> temp_idxA[iii] >> temp_idxB[iii];
717  // intersection_full_file.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
718  intersection_full_file >> temp_idxA[iii] >> temp_idxB[iii];
719  intersection_full_file.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
720  }
721  intersection_full_file.close();
722  }
723 
724  if(nodes > 1)
725  {
726  WorldComm.broadcast(nbOfIntersections);
727 
728  if(rank != 0)
729  {
730  temp_idxA.resize(nbOfIntersections);
731  temp_idxB.resize(nbOfIntersections);
732  }
733  WorldComm.broadcast(temp_idxA);
734  WorldComm.broadcast(temp_idxB);
735  }
736 
737  inter_mediator_A.reserve(nbOfIntersections);
738  inter_mediator_B.reserve(nbOfIntersections);
739 
740  int mediator_idx;
741  for(int iii = 0; iii < nbOfIntersections; ++iii)
742  {
743  mediator_idx = equivalence_table_system_to_mediator.at(temp_idxA[iii]);
744  inter_mediator_B.emplace(std::make_pair(mediator_idx, temp_idxB[iii]));
745  }
746 
747  int system_idx;
748  for(unsigned int iii = 0; iii < equivalence_table_system_to_mediator.size(); ++iii)
749  {
750  system_idx = equivalence_table_mediator_to_system.at(iii);
751  inter_mediator_A.emplace(std::make_pair(iii, system_idx));
752  }
753 };
void carl::set_intersection_tables ( const libMesh::Parallel::Communicator &  WorldComm,
const libMesh::Mesh &  mesh_intersection,
const std::string &  intersection_full_table_Filename,
const std::string &  equivalence_table_A_Filename,
const std::string &  equivalence_table_B_Filename,
const std::unordered_map< int, int > &  equivalence_table_A_to_R_A,
const std::unordered_map< int, int > &  equivalence_table_B_to_R_B,
std::unordered_map< int, std::pair< int, int > > &  full_intersection_pairs_map,
std::unordered_map< int, std::pair< int, int > > &  full_intersection_restricted_pairs_map,
std::unordered_map< int, int > &  local_intersection_meshI_to_inter_map 
)

Definition at line 595 of file mesh_tables.cpp.

609 {
610  // Start by reading and broadcasting the global intersection table
611  std::unordered_map<int,int> full_intersection_meshI_to_inter_map;
612  set_full_intersection_tables(WorldComm,intersection_full_table_Filename,
613  full_intersection_pairs_map,full_intersection_meshI_to_inter_map);
614 
615  // Set up the restricted intersection pairs table
616  set_restricted_intersection_pairs_table(full_intersection_pairs_map,
617  equivalence_table_A_to_R_A,equivalence_table_B_to_R_B,
618  full_intersection_restricted_pairs_map);
619 
620  // Build in each processor its own intersection table
621  libMesh::MeshBase::const_element_iterator elemIt = mesh_intersection.active_local_elements_begin();
622  const libMesh::MeshBase::const_element_iterator end_elemIt = mesh_intersection.active_local_elements_end();
623  int local_nbOfInters = mesh_intersection.n_active_local_elem();
624  local_intersection_meshI_to_inter_map.reserve(local_nbOfInters);
625 
626  // Some dummy indexes used
627  int idxI_table = -1;
628  int interID = -1;
629 
630  for( ; elemIt != end_elemIt; ++ elemIt)
631  {
632  const libMesh::Elem* elem = *elemIt;
633  idxI_table = elem->id();
634  interID = full_intersection_meshI_to_inter_map[idxI_table];
635  local_intersection_meshI_to_inter_map[idxI_table] = interID;
636  }
637 };
carl::set_restricted_intersection_pairs_table
void set_restricted_intersection_pairs_table(const std::unordered_map< int, std::pair< int, int > > &full_intersection_pairs_map, const std::unordered_map< int, int > &equivalence_table_A_to_R_A, const std::unordered_map< int, int > &equivalence_table_B_to_R_B, std::unordered_map< int, std::pair< int, int > > &full_intersection_restricted_pairs_map)
Definition: mesh_tables.cpp:402
carl::set_full_intersection_tables
void set_full_intersection_tables(const libMesh::Parallel::Communicator &WorldComm, const std::string &intersection_full_table_Filename, std::unordered_map< int, std::pair< int, int > > &full_intersection_pairs_map, std::unordered_map< int, int > &full_intersection_meshI_to_inter_map)
Definition: mesh_tables.cpp:437
void carl::set_local_intersection_tables ( const libMesh::Parallel::Communicator &  WorldComm,
const libMesh::Mesh &  mesh_intersection,
const std::string &  intersection_local_table_Filename,
const std::string &  equivalence_table_A_Filename,
const std::string &  equivalence_table_B_Filename,
const std::unordered_map< int, int > &  equivalence_table_A_to_R_A,
const std::unordered_map< int, int > &  equivalence_table_B_to_R_B,
std::unordered_map< int, std::pair< int, int > > &  local_intersection_pairs_map,
std::unordered_map< int, std::pair< int, int > > &  local_intersection_restricted_pairs_map,
std::unordered_map< int, int > &  local_intersection_meshI_to_inter_map 
)

Definition at line 639 of file mesh_tables.cpp.

653 {
654  // Start by reading and broadcasting the global intersection table
655  std::unordered_map<int,int> full_intersection_meshI_to_inter_map;
656  read_local_intersection_tables(WorldComm,intersection_local_table_Filename,
657  local_intersection_pairs_map,full_intersection_meshI_to_inter_map);
658 
659  // Set up the restricted intersection pairs table
660  set_restricted_intersection_pairs_table(local_intersection_pairs_map,
661  equivalence_table_A_to_R_A,equivalence_table_B_to_R_B,
662  local_intersection_restricted_pairs_map);
663 
664  // Build in each processor its own intersection table
665  libMesh::MeshBase::const_element_iterator elemIt = mesh_intersection.active_local_elements_begin();
666  const libMesh::MeshBase::const_element_iterator end_elemIt = mesh_intersection.active_local_elements_end();
667  int local_nbOfInters = mesh_intersection.n_active_local_elem();
668  local_intersection_meshI_to_inter_map.reserve(local_nbOfInters);
669 
670  // Some dummy indexes used
671  int idxI_table = -1;
672  int interID = -1;
673 
674  for( ; elemIt != end_elemIt; ++ elemIt)
675  {
676  const libMesh::Elem* elem = *elemIt;
677  idxI_table = elem->id();
678  interID = full_intersection_meshI_to_inter_map[idxI_table];
679  local_intersection_meshI_to_inter_map[idxI_table] = interID;
680  }
681 };
carl::set_restricted_intersection_pairs_table
void set_restricted_intersection_pairs_table(const std::unordered_map< int, std::pair< int, int > > &full_intersection_pairs_map, const std::unordered_map< int, int > &equivalence_table_A_to_R_A, const std::unordered_map< int, int > &equivalence_table_B_to_R_B, std::unordered_map< int, std::pair< int, int > > &full_intersection_restricted_pairs_map)
Definition: mesh_tables.cpp:402
carl::read_local_intersection_tables
void read_local_intersection_tables(const libMesh::Parallel::Communicator &WorldComm, const std::string &intersection_local_table_Filename, std::unordered_map< int, std::pair< int, int > > &local_intersection_pairs_map, std::unordered_map< int, int > &local_intersection_meshI_to_inter_map)
Definition: mesh_tables.cpp:546
void carl::set_restricted_intersection_pairs_table ( const std::unordered_map< int, std::pair< int, int > > &  full_intersection_pairs_map,
const std::unordered_map< int, int > &  equivalence_table_A_to_R_A,
const std::unordered_map< int, int > &  equivalence_table_B_to_R_B,
std::unordered_map< int, std::pair< int, int > > &  full_intersection_restricted_pairs_map 
)

Definition at line 402 of file mesh_tables.cpp.

407 {
408  // "Resize" the final output
409  intersection_restricted_pairs_map.reserve(equivalence_table_A_to_R_A.size());
410 
411  int interID = -1;
412  int idxA = -1;
413  int idxB = -1;
414 
415  int idxRA = -1;
416  int idxRB = -1;
417 
418  std::unordered_map<int,std::pair<int,int> >::const_iterator mapIt =
419  intersection_pairs_map.begin();
420  std::unordered_map<int,std::pair<int,int> >::const_iterator end_mapIt =
421  intersection_pairs_map.end();
422 
423  for( ; mapIt != end_mapIt; ++mapIt)
424  {
425  interID = mapIt->first;
426  idxA = mapIt->second.first;
427  idxB = mapIt->second.second;
428 
429  idxRA = equivalence_table_A_to_R_A.at(idxA);
430  idxRB = equivalence_table_B_to_R_B.at(idxB);
431 
432  intersection_restricted_pairs_map[interID] =
433  std::pair<int,int>(idxRA,idxRB);
434  }
435 }
void carl::set_weight_function_domain_idx ( std::string &  filename,
int &  domain_Idx_BIG,
int &  nb_of_domain_Idx_micro,
std::vector< int > &  domain_Idx_micro,
std::vector< int > &  domain_Idx_coupling 
)

Definition at line 4 of file mesh_tables.cpp.

10 {
11  std::ifstream dataF(filename);
12 
13  // Buffer string
14  std::string bufferLine;
15  std::stringstream dataBuffer;
16 
17  // Read info until the file ends
18  while(std::getline(dataF,bufferLine))
19  {
20  if(bufferLine.compare("$MacroDomainIdx")==0)
21  {
22  dataF >> domain_Idx_BIG;
23  }
24 
25  if(bufferLine.compare("$NbOfMicroDomainIdx")==0)
26  {
27  dataF >> nb_of_domain_Idx;
28 
29  domain_Idx_micro.resize(nb_of_domain_Idx);
30  domain_Idx_coupling.resize(nb_of_domain_Idx);
31  }
32 
33  if(bufferLine.compare("$MicroDomainIdxs")==0)
34  {
35  for(int iii = 0; iii < nb_of_domain_Idx; ++iii )
36  {
37  dataF >> domain_Idx_micro[iii];
38  }
39  }
40 
41  if(bufferLine.compare("$CouplingDomainIdxs")==0)
42  {
43  for(int iii = 0; iii < nb_of_domain_Idx; ++iii )
44  {
45  dataF >> domain_Idx_coupling[iii];
46  }
47  }
48  }
49 
50  dataF.close();
51 }
void carl::solve_linear_PETSC ( libMesh::PetscMatrix< libMesh::Number > &  A,
libMesh::PetscVector< libMesh::Number > &  b,
libMesh::PetscVector< libMesh::Number > &  x,
KSP &  ksp,
PC &  pc 
)

Definition at line 270 of file PETSC_matrix_operations.cpp.

274 {
275  /*
276  * Solve the system A*x = b using PETSc's linear solver, using a
277  * Krilov method, with the linear solver context "ksp" and
278  * preconditioner "pc".
279  */
280 
281  // Set up inital variables
282 }
int carl::voigt_index_converter ( int  aaa,
int  bbb 
)

Definition at line 93 of file common_functions.cpp.

94 {
95  if(aaa == bbb)
96  {
97  // 00 -> 0, 11 -> 1, 22 -> 2
98  return aaa;
99  }
100  else
101  {
102  // 12, 21 -> 3
103  // 02, 20 -> 4
104  // 01, 10 -> 5
105  int sum = aaa + bbb;
106 
107  if(sum == 3)
108  {
109  return 3;
110  }
111  else if(sum == 2)
112  {
113  return 4;
114  }
115  else if(sum == 1)
116  {
117  return 5;
118  }
119  }
120 
121  std::cerr << "Bad indexes! " << aaa << " " << bbb << std::endl;
122  homemade_error_msg(" You shouldn't be here! (voigt_index_converter)");
123  return -1;
124 };
homemade_error_msg
#define homemade_error_msg(msg)
Definition: common_header.h:73
void carl::write_PETSC_matrix ( Mat  input_mat,
const std::string &  filename,
int  rank,
MPI_Comm  comm = PETSC_COMM_WORLD,
int  dim = 1 
)

Definition at line 415 of file PETSC_matrix_operations.cpp.

417 {
418  PetscViewer viewer;
419  PetscViewerBinaryOpen(comm,filename.c_str(),FILE_MODE_WRITE,&viewer);
420  MatView(input_mat,viewer);
421 
422  PetscViewerDestroy(&viewer);
423 
424  // Hack to guarantee the proper matrix division without forcing libMesh's --enable-blocked-storage
425  if(rank == 0)
426  {
427  std::ofstream mat_info(filename + ".info");
428  mat_info << "-matload_block_size " << std::to_string(dim) << std::endl;
429  mat_info.close();
430  }
431 }
void carl::write_PETSC_matrix ( libMesh::PetscMatrix< libMesh::Number > &  input_mat,
const std::string &  filename,
int  dim = 1 
)

Definition at line 433 of file PETSC_matrix_operations.cpp.

435 {
436  write_PETSC_matrix(input_mat.mat(),filename,input_mat.comm().rank(),input_mat.comm().get(),dim);
437 }
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
void carl::write_PETSC_matrix_MATLAB ( Mat  input_mat,
const std::string &  filename,
MPI_Comm  comm = PETSC_COMM_WORLD 
)

Definition at line 356 of file PETSC_matrix_operations.cpp.

357 {
358  PetscViewer viewer;
359  PetscViewerCreate(comm,&viewer);
360  PetscViewerASCIIOpen(comm,filename.c_str(),&viewer);
361  PetscViewerPushFormat (viewer,PETSC_VIEWER_ASCII_MATLAB);
362  MatView(input_mat,viewer);
363 
364  PetscViewerDestroy(&viewer);
365 }
void carl::write_PETSC_vector ( libMesh::PetscVector< libMesh::Number > &  input_vec,
const std::string &  filename,
int  dim = 1 
)

Definition at line 338 of file PETSC_matrix_operations.cpp.

340 {
341  write_PETSC_vector(input_vec.vec(),filename,input_vec.comm().rank(),input_vec.comm().get(),dim);
342 }
carl::write_PETSC_vector
void write_PETSC_vector(libMesh::PetscVector< libMesh::Number > &input_vec, const std::string &filename, int dim=1)
Definition: PETSC_matrix_operations.cpp:338
void carl::write_PETSC_vector ( Vec  input_vec,
const std::string &  filename,
int  rank,
MPI_Comm  comm = PETSC_COMM_WORLD,
int  dim = 1 
)

Definition at line 320 of file PETSC_matrix_operations.cpp.

322 {
323  PetscViewer viewer;
324  PetscViewerBinaryOpen(comm,filename.c_str(),FILE_MODE_WRITE,&viewer);
325  VecView(input_vec,viewer);
326 
327  PetscViewerDestroy(&viewer);
328 
329  // Hack to guarantee the proper vector division without forcing libMesh's --enable-blocked-storage
330  if(rank == 0)
331  {
332  std::ofstream vec_info(filename + ".info");
333  vec_info << "-vecload_block_size " << std::to_string(dim) << std::endl;
334  vec_info.close();
335  }
336 }
void carl::write_PETSC_vector_MATLAB ( Vec  input_vec,
const std::string &  filename,
MPI_Comm  comm = PETSC_COMM_WORLD 
)

Definition at line 344 of file PETSC_matrix_operations.cpp.

346 {
347  PetscViewer viewer;
348  PetscViewerCreate(comm,&viewer);
349  PetscViewerASCIIOpen(comm,filename.c_str(),&viewer);
350  PetscViewerPushFormat (viewer,PETSC_VIEWER_ASCII_MATLAB);
351  VecView(input_vec,viewer);
352 
353  PetscViewerDestroy(&viewer);
354 }