Add time scheme for mechanics

.clang-format

@@ -0,0 +1,47 @@
+BasedOnStyle: LLVM
+---
+Language: Cpp
+Standard: Cpp11
+IndentWidth: 4
+TabWidth: 4
+UseTab: Never
+ColumnLimit: 100
+AlwaysBreakTemplateDeclarations: Yes
+BreakConstructorInitializers : BeforeColon
+ConstructorInitializerAllOnOneLineOrOnePerLine: true
+FixNamespaceComments : true
+IncludeBlocks: Regroup
+IncludeCategories:
+  # aster main header
+  - Regex:           '^"aster[cx]*.h"$'
+    Priority:        -1
+  # Headers in "" with extension and without /.
+  - Regex:           '"([a-z0-9\._])+"'
+    Priority:        3
+  # Headers in "" with extension and with /.
+  - Regex:           '"([a-z0-9\./\-_])+"'
+    Priority:        4
+  # Headers in <> without extension.
+  - Regex:           '<([a-z0-9/\-_])+>'
+    Priority:        5
+  # Headers in <> with extension.
+  - Regex:           '<([a-z0-9\./\-_])+>'
+    Priority:        6
+
+# not properly applied by clang-format-11, will be set 'true' when updated
+Cpp11BracedListStyle : false
+SpaceBeforeCpp11BracedList: true
+
+SpacesInAngles: true
+SpacesInSquareBrackets: false
+
+SpaceAfterCStyleCast: false
+SpacesInParentheses: true
+SpacesInCStyleCastParentheses: false
+
+# for clang-format>=17
+# SpacesInParens: Custom
+# SpacesInParensOptions:
+#   InConditionalStatements: true
+#   InEmptyParentheses: false
+#   InCStyleCasts: false

apps/acoustic_eigs/src/acoustic_eigs.cpp

@@ -9,14 +9,13 @@
  */
 
 #include <cstddef>
+#include <filesystem>
+#include <iomanip>
 #include <iostream>
 #include <regex>
 #include <set>
 #include <string>
-#include <filesystem>
-#include <iomanip>
-
-
+#include <unistd.h>
 
 #include "sol/sol.hpp"
 #include "diskpp/common/eigen.hpp"
@@ -72,7 +71,7 @@ auto
 acoustic_eigs_dg(Mesh& msh, size_t degree,
     const typename Mesh::coordinate_type eta,
     disk::silo_database& silo)
-{   
+{
     std::cout << "DG eigsolver" << std::endl;
     auto cvf = connectivity_via_faces(msh);
     using T = typename Mesh::coordinate_type;
@@ -83,13 +82,13 @@ acoustic_eigs_dg(Mesh& msh, size_t degree,
 
     auto cbs = disk::scalar_basis_size(degree, Mesh::dimension);
     auto assm = make_discontinuous_galerkin_eigenvalue_assembler(msh, cbs);
-    
+
     timecounter tc;
     tc.tic();
     for (auto& tcl : msh)
     {
         auto tbasis = disk::basis::scaled_monomial_basis(msh, tcl, degree, basis_rescaling);
-        
+
         matrix_type M = integrate(msh, tcl, tbasis, tbasis);
         matrix_type K = integrate(msh, tcl, grad(tbasis), grad(tbasis));
 
@@ -98,15 +97,15 @@ acoustic_eigs_dg(Mesh& msh, size_t degree,
 
         auto fcs = faces(msh, tcl);
         for (auto& fc : fcs)
-        {   
+        {
             auto n     = normal(msh, tcl, fc);
             auto eta_l = eta / diameter(msh, fc);
-            
+
             auto nv = cvf.neighbour_via(msh, tcl, fc);
             if (nv) {
                 matrix_type Att = matrix_type::Zero(tbasis.size(), tbasis.size());
                 matrix_type Atn = matrix_type::Zero(tbasis.size(), tbasis.size());
-                
+
                 auto ncl = nv.value();
                 auto nbasis = disk::basis::scaled_monomial_basis(msh, ncl, degree, basis_rescaling);
                 assert(tbasis.size() == nbasis.size());
@@ -128,7 +127,7 @@ acoustic_eigs_dg(Mesh& msh, size_t degree,
                 //Att += - integrate(msh, fc, grad(tbasis).dot(n), tbasis);
                 //Att += - integrate(msh, fc, tbasis, grad(tbasis).dot(n));
                 //assm.assemble(msh, tcl, tcl, Att);
-            }   
+            }
         }
     }
 
@@ -183,7 +182,7 @@ acoustic_eigs_dg(Mesh& msh, size_t degree,
             auto ofs = cbs * offset(msh, cl);
             u.push_back(eigvecs(ofs, col));
         }
-        
+
         std::string vname = "eigfun_" + std::to_string(col);
         silo.add_variable("mesh", vname, u, disk::zonal_variable_t);
     }
@@ -194,7 +193,7 @@ void solve_feast_dense(auto assm, disk::dynamic_matrix<T>& eigvecs,
     disk::dynamic_vector<T>& eigvals)
 {
     timecounter tc;
-
+#ifdef HAVE_MUMPS
     std::cout << "MUMPS factorization..." << std::flush;
     tc.tic();
     disk::solvers::mumps_solver<T> AFF_lu_mumps;
@@ -218,6 +217,9 @@ void solve_feast_dense(auto assm, disk::dynamic_matrix<T>& eigvecs,
     disk::solvers::feast(params, KTT, assm.BTT, eigvecs, eigvals);
 
     std::cout << "Eigensolver time: " << tc.toc() << " seconds\n";
+#else
+    std::cerr << "MUMPS is needed" << std::endl;
+#endif
 }
 
 template<typename T>
@@ -226,6 +228,7 @@ void solve_feast_mf(auto assm, disk::dynamic_matrix<T>& eigvecs,
 {
     timecounter tc;
 
+#ifdef HAVE_PARDISO
     Eigen::PardisoLDLT< Eigen::SparseMatrix<T> > AFF_lu(assm.AFF);
 
     auto apply_A = [&]<int ncols>(
@@ -234,7 +237,6 @@ void solve_feast_mf(auto assm, disk::dynamic_matrix<T>& eigvecs,
         Eigen::Matrix<T, Eigen::Dynamic, ncols> z = assm.AFT*v;
         return assm.ATT*v - assm.ATF*AFF_lu.solve(z);
     };
-
     std::cout << "FEAST eigensolver (matrix-free)" << std::endl;
     tc.tic();
     disk::solvers::feast_eigensolver_params<T> params;
@@ -247,6 +249,9 @@ void solve_feast_mf(auto assm, disk::dynamic_matrix<T>& eigvecs,
     disk::solvers::feast_mf(params, apply_A, assm.BTT, eigvecs, eigvals);
 
     std::cout << "Eigensolver time: " << tc.toc() << " seconds\n";
+#else
+    std::cerr << "Pardiso is needed" << std::endl;
+#endif
 }
 
 
@@ -256,6 +261,7 @@ void solve_bjd_mf(auto assm, disk::dynamic_matrix<T>& eigvecs,
 {
     timecounter tc;
 
+#ifdef HAVE_MUMPS
     //Eigen::PardisoLDLT< Eigen::SparseMatrix<T> > AFF_lu(assm.AFF);
     disk::solvers::mumps_solver<T> AFF_lu;
     AFF_lu.symmetric(true);
@@ -280,6 +286,9 @@ void solve_bjd_mf(auto assm, disk::dynamic_matrix<T>& eigvecs,
     disk::solvers::block_jacobi_davidson(params, apply_A,
         assm.BTT, eigvecs, eigvals);
     std::cout << "Eigensolver time: " << tc.toc() << " seconds\n";
+#else
+    std::cerr << "MUMPS is needed" << std::endl;
+#endif
 }
 
 #if 0
@@ -397,7 +406,7 @@ acoustic_eigs_hho(const Mesh& msh, const config& cfg, disk::silo_database& silo)
             auto ofs = cbasis_type::size_of_degree(di.cell) * offset(msh, cl);
             u.push_back(eigvecs(ofs, col));
         }
-        
+
         std::string vname = "hho_eigfun_" + std::to_string(col);
         silo.add_variable("hmesh", vname, u, disk::zonal_variable_t);
     }
@@ -442,7 +451,7 @@ int main(int argc, char **argv)
         case 'f':
             cfg.mesh_filename = optarg;
             break;
-        
+
         case 'k':
             cfg.order = std::stoul(optarg);
             break;
@@ -475,30 +484,37 @@ int main(int argc, char **argv)
 
 
     if (cfg.mesh_filename != "") {
- 
+
         if (std::regex_match(cfg.mesh_filename, std::regex(".*\\.geo2s$") ))
         {
             std::cout << "Guessed mesh format: GMSH 2D simplicials" << std::endl;
             using mesh_type = disk::triangular_mesh<T>;
             mesh_type msh;
+#ifdef HAVE_GMSH
             disk::gmsh_geometry_loader< mesh_type > loader;
             loader.read_mesh(cfg.mesh_filename);
             loader.populate_mesh(msh);
-
             run_eigsolver(msh, cfg);
+#else
+    std::cerr << "GMSH is needed" << std::endl;
+#endif
             return 0;
         }
-  
+
         if (std::regex_match(cfg.mesh_filename, std::regex(".*\\.geo3s$") ))
         {
             std::cout << "Guessed mesh format: GMSH 3D simplicials" << std::endl;
             using mesh_type = disk::tetrahedral_mesh<T>;
             mesh_type msh;
+#ifdef HAVE_GMSH
             disk::gmsh_geometry_loader< mesh_type > loader;
             loader.read_mesh(cfg.mesh_filename);
             loader.populate_mesh(msh);
 
             run_eigsolver(msh, cfg);
+#else
+    std::cerr << "GMSH is needed" << std::endl;
+#endif
             return 0;
         }
     }
@@ -515,7 +531,7 @@ int main(int argc, char **argv)
 
         for (int i = 0; i < cfg.reflevels; i++) {
             mesher.refine();
-        
+
             std::cout << ">>>>>>>> DIAM: " << disk::average_diameter(msh) << std::endl;
             run_eigsolver(msh, cfg);
         }
@@ -533,7 +549,7 @@ int main(int argc, char **argv)
 
         for (int i = 0; i < cfg.reflevels; i++) {
             mesher.refine();
-        
+
             std::cout << ">>>>>>>> DIAM: " << disk::average_diameter(msh) << std::endl;
             run_eigsolver(msh, cfg);
         }
@@ -551,7 +567,7 @@ int main(int argc, char **argv)
 
         for (int i = 0; i < cfg.reflevels; i++) {
             mesher.refine();
-        
+
             std::cout << ">>>>>>>> DIAM: " << disk::average_diameter(msh) << std::endl;
             run_eigsolver(msh, cfg);
         }
@@ -570,7 +586,7 @@ int main(int argc, char **argv)
 
         for (int i = 0; i < cfg.reflevels; i++) {
             mesher.refine();
-        
+
             std::cout << ">>>>>>>> DIAM: " << disk::average_diameter(msh) << std::endl;
             run_eigsolver(msh, cfg);
         }
@@ -589,11 +605,11 @@ int main(int argc, char **argv)
             msh.transform( [&](const typename mesh_type::point_type& pt) {
                 return typename mesh_type::point_type{pt.x(), 1.1*pt.y()};
             } );
-        
+
             std::cout << ">>>>>>>> DIAM: " << disk::average_diameter(msh) << std::endl;
             run_eigsolver(msh, cfg);
         }
     }
-    
+
     return 0;
 }

apps/contact/src/common.hpp

@@ -1903,7 +1903,7 @@ class diffusion_condensed_assembler_nitsche_cells
             {
                 auto fb = make_scalar_monomial_basis(msh, fc, di.face_degree());
                 auto dirichlet_bf = m_bnd.dirichlet_boundary_func(face_id);
-                Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb, di.face_degree());
+                Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb);
                 Matrix<T, Dynamic, 1> rhs = make_rhs(msh, fc, fb, dirichlet_bf, di.face_degree());
                 //ret.block(face_i*fbs, 0, fbs, 1) = mass.llt().solve(rhs);
             }
@@ -2249,7 +2249,7 @@ class diffusion_full_assembler
             if (m_bnd.is_dirichlet_face( face_id))
             {
                 auto fb = disk::make_scalar_monomial_basis(msh, fc, di.face_degree());
-                Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb, di.face_degree());
+                Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb);
                 auto velocity = m_bnd.dirichlet_boundary_func(face_id);
                 Matrix<T, Dynamic, 1> rhs = make_rhs(msh, fc, fb, velocity, di.face_degree());
                 svel.block(cbs + i * fbs, 0, fbs, 1) = mass.llt().solve(rhs);
@@ -2569,7 +2569,7 @@ class diffusion_mix_full_assembler
             if (m_bnd.is_dirichlet_face( face_id))
             {
                 auto fb = disk::make_scalar_monomial_basis(msh, fc, di.face_degree());
-                Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb, di.face_degree());
+                Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb);
                 auto velocity = m_bnd.dirichlet_boundary_func(face_id);
                 Matrix<T, Dynamic, 1> rhs = make_rhs(msh, fc, fb, velocity, di.face_degree());
                 svel.block(cbs + i * fbs, 0, fbs, 1) = mass.llt().solve(rhs);
@@ -2893,7 +2893,7 @@ class contact_full_assembler
             if (dirichlet)
             {
                 auto fb = disk::make_scalar_monomial_basis(msh, fc, di.face_degree());
-                Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb, di.face_degree());
+                Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, fc, fb);
                 auto velocity = m_bnd.dirichlet_boundary_func(face_id);
                 Matrix<T, Dynamic, 1> rhs = make_rhs(msh, fc, fb, velocity);//, di.face_degree());
                 svel.block(cbs + i * fbs, 0, fbs, 1) = mass.llt().solve(rhs);
@@ -3119,7 +3119,7 @@ class contact_full_assembler_new
             auto fb = make_scalar_monomial_basis(msh, face, di.face_degree());
             auto dirichlet_fun  = m_bnd.dirichlet_boundary_func(face_id);
 
-            matrix_type mass = make_mass_matrix(msh, face, fb);// di.face_degree());
+            matrix_type mass = make_mass_matrix(msh, face, fb);
             vector_type rhs  = make_rhs(msh, face, fb, dirichlet_fun);// di.face_degree());
 
             sol.block(face_ofs,  0, fbs, 1) = mass.llt().solve(rhs);
@@ -3448,7 +3448,7 @@ class contact_face_assembler_new
                         Matrix<T, Dynamic, 1> robin = make_rhs(msh,bfc,fb,bnd.robin_boundary_func(face_id), face_degree);
                         assert (robin.size() == num_face_dofs);
 
-                        Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, bfc, fb, face_degree);
+                        Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, bfc, fb);
 
                         for (size_t i = 0; i < num_face_dofs; i++)
                         {
@@ -3518,7 +3518,7 @@ class contact_face_assembler_new
             auto fb = make_scalar_monomial_basis(msh, face, di.face_degree());
             auto dirichlet_fun  = m_bnd.dirichlet_boundary_func(face_id);
 
-            Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, face, fb, di.face_degree());
+            Matrix<T, Dynamic, Dynamic> mass = make_mass_matrix(msh, face, fb);
             Matrix<T, Dynamic, 1> rhs  = make_rhs(msh, face, fb, dirichlet_fun, di.face_degree());
 
             sol.block(face_ofs,  0, fb.size(), 1) = mass.llt().solve(rhs);

apps/contact/src/diffusion_nitsche_solver.hpp

@@ -206,7 +206,7 @@ run_hho_diffusion_nitsche_faces(const Mesh& msh,
         auto diff = realsol - fullsol;
         H1_error += diff.dot(A*diff);
 
-        matrix_type mass  = make_mass_matrix(msh, cl, cb, hdi.cell_degree());
+        matrix_type mass  = make_mass_matrix(msh, cl, cb);
         vector_type u_diff = diff.block(0, 0, cbs, 1);
         L2_error += u_diff.dot(mass * u_diff);
 
@@ -370,7 +370,7 @@ run_hho_diffusion_nitsche_cells_full(const Mesh& msh,
         auto diff = realsol - fullsol;
         H1_error += diff.dot(A*diff);
 
-        matrix_type mass  = make_mass_matrix(msh, cl, cb, hdi.cell_degree());
+        matrix_type mass  = make_mass_matrix(msh, cl, cb);
 
         vector_type u_diff = diff.block(0, 0, cbs, 1);
         L2_error += u_diff.dot(mass * u_diff);