diff --git a/test/test_model/CMakeLists.txt b/test/test_model/CMakeLists.txt index 2d7e0230f..a0c8081d7 100644 --- a/test/test_model/CMakeLists.txt +++ b/test/test_model/CMakeLists.txt @@ -1,30 +1,29 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # # @date creation: Fri Sep 03 2010 # @date last modification: Tue Jan 30 2018 # # @brief configuration for model tests # # @section LICENSE # # Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_subdirectory(patch_tests) -add_akantu_test(test_model_solver "Test for the solvers") +add_akantu_test(test_common "Tests the common part of the models") add_akantu_test(test_solid_mechanics_model "Test for the solid mechanics model") add_akantu_test(test_structural_mechanics_model "Test for the structural mechanics model") -add_akantu_test(test_non_local_toolbox "Test of the functionalities in the non-local toolbox") add_akantu_test(test_heat_transfer_model "Test heat transfer model") diff --git a/test/test_model/test_common/CMakeLists.txt b/test/test_model/test_common/CMakeLists.txt new file mode 100644 index 000000000..5d44f51bc --- /dev/null +++ b/test/test_model/test_common/CMakeLists.txt @@ -0,0 +1,2 @@ +add_akantu_test(test_model_solver "Test for the solvers") +add_akantu_test(test_non_local_toolbox "Test of the functionalities in the non-local toolbox") diff --git a/test/test_model/test_model_solver/CMakeLists.txt b/test/test_model/test_common/test_model_solver/CMakeLists.txt similarity index 100% rename from test/test_model/test_model_solver/CMakeLists.txt rename to test/test_model/test_common/test_model_solver/CMakeLists.txt diff --git a/test/test_model/test_model_solver/test_dof_manager_default.cc b/test/test_model/test_common/test_model_solver/test_dof_manager_default.cc similarity index 100% rename from test/test_model/test_model_solver/test_dof_manager_default.cc rename to test/test_model/test_common/test_model_solver/test_dof_manager_default.cc diff --git a/test/test_model/test_model_solver/test_dof_manager_default.verified b/test/test_model/test_common/test_model_solver/test_dof_manager_default.verified similarity index 100% rename from test/test_model/test_model_solver/test_dof_manager_default.verified rename to test/test_model/test_common/test_model_solver/test_dof_manager_default.verified diff --git a/test/test_model/test_model_solver/test_model_solver.cc b/test/test_model/test_common/test_model_solver/test_model_solver.cc similarity index 100% rename from test/test_model/test_model_solver/test_model_solver.cc rename to test/test_model/test_common/test_model_solver/test_model_solver.cc diff --git a/test/test_model/test_model_solver/test_model_solver.py b/test/test_model/test_common/test_model_solver/test_model_solver.py similarity index 100% rename from test/test_model/test_model_solver/test_model_solver.py rename to test/test_model/test_common/test_model_solver/test_model_solver.py diff --git a/test/test_model/test_model_solver/test_model_solver.verified b/test/test_model/test_common/test_model_solver/test_model_solver.verified similarity index 100% rename from test/test_model/test_model_solver/test_model_solver.verified rename to test/test_model/test_common/test_model_solver/test_model_solver.verified diff --git a/test/test_model/test_model_solver/test_model_solver_dynamic.cc b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc similarity index 100% rename from test/test_model/test_model_solver/test_model_solver_dynamic.cc rename to test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc diff --git a/test/test_model/test_model_solver/test_model_solver_dynamic.py b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.py similarity index 100% rename from test/test_model/test_model_solver/test_model_solver_dynamic.py rename to test/test_model/test_common/test_model_solver/test_model_solver_dynamic.py diff --git a/test/test_model/test_model_solver/test_model_solver_dynamic_explicit.verified b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic_explicit.verified similarity index 100% rename from test/test_model/test_model_solver/test_model_solver_dynamic_explicit.verified rename to test/test_model/test_common/test_model_solver/test_model_solver_dynamic_explicit.verified diff --git a/test/test_model/test_model_solver/test_model_solver_dynamic_implicit.verified b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic_implicit.verified similarity index 100% rename from test/test_model/test_model_solver/test_model_solver_dynamic_implicit.verified rename to test/test_model/test_common/test_model_solver/test_model_solver_dynamic_implicit.verified diff --git a/test/test_model/test_model_solver/test_model_solver_dynamic_petsc.cc b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic_petsc.cc similarity index 98% rename from test/test_model/test_model_solver/test_model_solver_dynamic_petsc.cc rename to test/test_model/test_common/test_model_solver/test_model_solver_dynamic_petsc.cc index 41f3b40c4..eba74e8ea 100644 --- a/test/test_model/test_model_solver/test_model_solver_dynamic_petsc.cc +++ b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic_petsc.cc @@ -1,838 +1,838 @@ /** * @file test_model_solver_dynamic.cc * * @author Nicolas Richart * * @date creation: Wed Apr 13 2016 * @date last modification: Tue Feb 20 2018 * * @brief Test default dof manager * * @section LICENSE * * Copyright (©) 2016-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "communicator.hh" #include "element_group.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "non_linear_solver.hh" /* -------------------------------------------------------------------------- */ #include "boundary_condition_functor.hh" #include "mpi_communicator_data.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #include "dumper_element_partition.hh" #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #ifndef EXPLICIT #define EXPLICIT true #endif template void CHECK_ERR_CXX(func && func_, PetscErrorCode ierr) { if (PetscUnlikely(ierr != 0)) { const char * desc; PetscErrorMessage(ierr, &desc, nullptr); AKANTU_EXCEPTION("Error in PETSc call to \'" << func_ << "\': " << desc); } } using namespace akantu; static void genMesh(Mesh & mesh, UInt nb_nodes); class MyModel { public: MyModel(Real F, Mesh & mesh, bool lumped) : nb_dofs(mesh.getNbNodes()), nb_elements(mesh.getNbElement(_segment_2)), lumped(lumped), E(1.), A(1.), rho(1.), mesh(mesh), displacement(nb_dofs, 1, "disp"), velocity(nb_dofs, 1, "velo"), acceleration(nb_dofs, 1, "accel"), blocked(nb_dofs, 1, "blocked"), forces(nb_dofs, 1, "force_ext"), internal_forces(nb_dofs, 1, "force_int"), stresses(nb_elements, 1, "stress"), strains(nb_elements, 1, "strain"), initial_lengths(nb_elements, 1, "L0") { auto n_global = mesh.getNbGlobalNodes(); int n_local = 0; std::vector nodes_global_ids(nb_dofs); for (auto && data : enumerate(nodes_global_ids)) { auto n = std::get<0>(data); n_local += mesh.isLocalOrMasterNode(n); std::get<1>(data) = mesh.getNodeGlobalId(n); } mpi_comm = dynamic_cast( mesh.getCommunicator().getCommunicatorData()) .getMPICommunicator(); MeshAccessor mesh_accessor(mesh); ierr = ISLocalToGlobalMappingCreate( mpi_comm, 1, mesh.getNbNodes(), nodes_global_ids.data(), PETSC_COPY_VALUES, &petsc_local_to_global); CHECK_ERR_CXX("ISLocalToGlobalMappingCreate", ierr); auto setName = [](auto && Obj, auto && name) { PetscObjectSetName(reinterpret_cast(Obj), name); }; ierr = VecCreate(mpi_comm, &rhs); ierr = VecSetSizes(rhs, n_local, n_global); ierr = VecSetFromOptions(rhs); ierr = VecSetLocalToGlobalMapping(rhs, petsc_local_to_global); setName(rhs, "rhs"); ierr = VecDuplicate(rhs, &x); ierr = VecDuplicate(rhs, &x_save); ierr = VecDuplicate(rhs, &dx); ierr = VecDuplicate(rhs, &f_int); ierr = VecDuplicate(rhs, &f_dirichlet); setName(x, "x"); setName(x_save, "x save"); setName(dx, "dx"); setName(f_int, "f_int"); setName(f_dirichlet, "f_dirichlet"); ierr = MatCreate(mpi_comm, &M); ierr = MatSetSizes(M, n_local, n_local, n_global, n_global); ierr = MatSetFromOptions(M); ierr = MatSetOption(M, MAT_SYMMETRIC, PETSC_TRUE); ierr = MatSetOption(M, MAT_ROW_ORIENTED, PETSC_TRUE); ierr = MatSetUp(M); ierr = MatSetLocalToGlobalMapping(M, petsc_local_to_global, petsc_local_to_global); setName(M, "M"); assembleMass(); ierr = MatDuplicate(M, MAT_DO_NOT_COPY_VALUES, &K); setName(K, "K"); ierr = MatDuplicate(M, MAT_DO_NOT_COPY_VALUES, &J); setName(J, "J"); ierr = SNESCreate(mpi_comm, &snes); ierr = SNESSetFromOptions(snes); ierr = SNESSetFunction(snes, rhs, MyModel::FormFunction, this); ierr = SNESSetJacobian(snes, J, J, MyModel::FormJacobian, this); PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INDEX); displacement.set(0.); velocity.set(0.); acceleration.set(0.); forces.set(0.); blocked.set(false); blocked(0, 0) = true; blocked(nb_dofs - 1, 0) = true; displacement(0, 0) = 0; displacement(nb_dofs - 1, 0) = 1; for (auto && data : zip(make_view(this->mesh.getConnectivity(_segment_2), 2), make_view(this->initial_lengths))) { const auto & conn = std::get<0>(data); auto & L = std::get<1>(data); auto p1 = this->mesh.getNodes()(conn(0), _x); auto p2 = this->mesh.getNodes()(conn(1), _x); L = std::abs(p2 - p1); } } // static PetscErrorCode SNESMonitor(SNES snes,PetscInt its,PetscReal // fnorm,void *ctx) { // auto & _this = *reinterpret_cast(ctx); // //SNESMonitorDefault(snes, its, fnorm, PETSC_VIEWER_STDOUT_WORLD); // } static PetscErrorCode FormFunction(SNES /*snes*/, Vec /*dx*/, Vec /*f*/, void * ctx) { auto & _this = *reinterpret_cast(ctx); _this.assembleResidual(); return 0; } static PetscErrorCode FormJacobian(SNES /*snes*/, Vec /*dx*/, Mat /*J*/, Mat /*P*/, void * ctx) { auto & _this = *reinterpret_cast(ctx); _this.assembleJacobian(); return 0; } ~MyModel() { ierr = MatDestroy(&M); ierr = MatDestroy(&K); ierr = MatDestroy(&J); ierr = VecDestroy(&rhs); ierr = VecDestroy(&x); ierr = VecDestroy(&dx); ierr = VecDestroy(&x_save); ierr = VecDestroy(&f_int); PetscFinalize(); } void solveStep() { std::cout << "solveStep" << std::endl; copy(x_save, displacement); ierr = SNESSolve(snes, NULL, dx); CHECK_ERR_CXX("SNESSolve", ierr); setSolutionToDisplacement(); assembleResidual(); } void applyBC() { std::vector rows; for (auto && data : enumerate(blocked)) { if (std::get<1>(data)) { rows.push_back(std::get<0>(data)); } } copy(x, displacement); ierr = MatZeroRowsColumnsLocal(J, rows.size(), rows.data(), 1., x, f_dirichlet); VecView(f_dirichlet, PETSC_VIEWER_STDOUT_WORLD); CHECK_ERR_CXX("MatZeroRowsColumnsLocal", ierr); } void setSolutionToDisplacement() { std::cout << "setSolutionToDisplacement" << std::endl; ierr = VecWAXPY(x, 1, x_save, dx); copy(displacement, x); } void assembleJacobian() { std::cout << "assembleJacobian" << std::endl; setSolutionToDisplacement(); assembleStiffness(); ierr = MatZeroEntries(J); CHECK_ERR_CXX("MatZeroEntries", ierr); ierr = MatAXPY(J, 1., K, SAME_NONZERO_PATTERN); CHECK_ERR_CXX("MatAXPY", ierr); MatView(J, PETSC_VIEWER_STDOUT_WORLD); applyBC(); MatView(J, PETSC_VIEWER_STDOUT_WORLD); } void assembleMass() { std::cout << "assembleMass" << std::endl; ierr = MatZeroEntries(M); CHECK_ERR_CXX("MatZeroEntries", ierr); Array m_all_el(this->nb_elements, 4); Matrix m(2, 2); m(0, 0) = m(1, 1) = 2; m(0, 1) = m(1, 0) = 1; // under integrated // m(0, 0) = m(1, 1) = 3./2.; // m(0, 1) = m(1, 0) = 3./2.; // lumping the mass matrix // m(0, 0) += m(0, 1); // m(1, 1) += m(1, 0); // m(0, 1) = m(1, 0) = 0; for (auto && data : zip(make_view(this->mesh.getConnectivity(_segment_2), 2), make_view(m_all_el, 2, 2))) { const auto & conn = std::get<0>(data); auto & m_el = std::get<1>(data); UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); m_el = m; m_el *= rho * A * L / 6.; Vector conn_int(conn.size()); for (auto && data : zip(conn_int, conn)) { std::get<0>(data) = std::get<1>(data); } ierr = MatSetValuesLocal(M, conn_int.size(), conn_int.storage(), conn_int.size(), conn_int.storage(), m.storage(), ADD_VALUES); } ierr = MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY); ierr = MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY); ierr = MatSetOption(M, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE); PetscViewer viewer; ierr = PetscViewerASCIIOpen(mpi_comm, "M.mtx", &viewer); PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATRIXMARKET); ierr = MatView(M, viewer); PetscViewerPopFormat(viewer); ierr = PetscViewerDestroy(&viewer); // this->getDOFManager().assembleElementalMatricesToMatrix( // "M", "disp", m_all_el, _segment_2); is_mass_assembled = true; } // MatrixType getMatrixType(const ID &) { return _symmetric; } // void assembleMatrix(const ID & matrix_id) { // if (matrix_id == "K") { // if (not is_stiffness_assembled) // this->assembleStiffness(); // } else if (matrix_id == "M") { // if (not is_mass_assembled) // this->assembleMass(); // } else if (matrix_id == "C") { // // pass, no damping matrix // } else { // AKANTU_EXCEPTION("This solver does not know what to do with a matrix " // << matrix_id); // } // } void assembleLumpedMatrix(const ID & matrix_id) { std::cout << "assembleLumpedMatrix" << std::endl; AKANTU_EXCEPTION("This solver does not know what to do with a matrix " << matrix_id); } void assembleStiffness() { std::cout << "assembleStiffness" << std::endl; // SparseMatrix & K = this->getDOFManager().getMatrix("K"); // K.clear(); ierr = MatZeroEntries(K); CHECK_ERR_CXX("MatZeroEntries", ierr); Matrix k(2, 2); k(0, 0) = k(1, 1) = 1; k(0, 1) = k(1, 0) = -1; Array k_all_el(this->nb_elements, 4); auto k_it = k_all_el.begin(2, 2); auto cit = this->mesh.getConnectivity(_segment_2).begin(2); auto cend = this->mesh.getConnectivity(_segment_2).end(2); for (; cit != cend; ++cit, ++k_it) { const auto & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); auto & k_el = *k_it; k_el = k; k_el *= E * A / L; Vector conn_int(conn.size()); for (auto && data : zip(conn_int, conn)) { std::get<0>(data) = std::get<1>(data); } ierr = MatSetValuesLocal(K, conn_int.size(), conn_int.storage(), conn_int.size(), conn_int.storage(), k_el.storage(), ADD_VALUES); } ierr = MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY); CHECK_ERR_CXX("MatAssemblyBegin", ierr); ierr = MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY); CHECK_ERR_CXX("MatAssemblyEnd", ierr); ierr = MatSetOption(K, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE); CHECK_ERR_CXX("MatSetOption", ierr); PetscViewer viewer; ierr = PetscViewerASCIIOpen(mpi_comm, "K.mtx", &viewer); CHECK_ERR_CXX("PetscViewerASCIIOpen", ierr); PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATRIXMARKET); ierr = MatView(K, viewer); CHECK_ERR_CXX("MatView", ierr); PetscViewerPopFormat(viewer); ierr = PetscViewerDestroy(&viewer); CHECK_ERR_CXX("PetscViewerDestroy", ierr); // this->getDOFManager().assembleElementalMatricesToMatrix( // "K", "disp", k_all_el, _segment_2); is_stiffness_assembled = true; } void copy(Array & y, Vec x) { std::cout << "copy <-" << std::endl; const PetscScalar * x_local; ierr = VecGetArrayRead(x, &x_local); for (auto && data : zip(y, range(x_local + 0, x_local + y.size()))) { std::get<0>(data) = std::get<1>(data); } ierr = VecRestoreArrayRead(x, &x_local); // VecView(x, PETSC_VIEWER_STDOUT_WORLD); // std::cout << y.getID() << " " << Vector(y.storage(), y.size()) // << std::endl; } void print(const Array & x) const { std::cout << x.getID() << " " << Vector(x.storage(), x.size()) << std::endl; } void copy(Vec x, const Array & y) { std::cout << "copy ->" << std::endl; PetscScalar * x_local; ierr = VecGetArray(x, &x_local); for (auto && data : zip(y, range(x_local + 0, x_local + y.size()))) { std::get<1>(data) = std::get<0>(data); } ierr = VecRestoreArray(x, &x_local); // std::cout << y.getID() << " " << Vector(y.storage(), y.size()) // << std::endl; // VecView(x, PETSC_VIEWER_STDOUT_WORLD); } void assembleResidual() { std::cout << "assembleResidual" << std::endl; // this->getDOFManager().assembleToResidual("disp", forces); setSolutionToDisplacement(); copy(rhs, forces); // VecAXPY(rhs, -1., f_dirichlet); print(displacement); this->assembleResidual(_not_ghost); // this->synchronize(SynchronizationTag::_user_1); // this->getDOFManager().assembleToResidual("disp", internal_forces, -1.); VecAXPY(rhs, 1., f_int); for (auto && data : enumerate(blocked)) { if(std::get<1>(data)) { VecSetValueLocal(rhs, std::get<0>(data), 0., INSERT_VALUES); } } VecAssemblyBegin(rhs); VecAssemblyEnd(rhs); VecView(rhs, PETSC_VIEWER_STDOUT_WORLD); } void assembleResidual(const GhostType & ghost_type) { std::cout << "assembleResidual" << std::endl; VecZeroEntries(f_int); auto cit = this->mesh.getConnectivity(_segment_2, ghost_type).begin(2); auto cend = this->mesh.getConnectivity(_segment_2, ghost_type).end(2); auto strain_it = this->strains.begin(); auto stress_it = this->stresses.begin(); auto L_it = this->initial_lengths.begin(); for (; cit != cend; ++cit, ++strain_it, ++stress_it, ++L_it) { const auto & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real u1 = this->displacement(n1, _x); Real u2 = this->displacement(n2, _x); *strain_it = (u2 - u1) / *L_it; *stress_it = E * *strain_it; Real f_n = A * *stress_it; std::cout << n1 << "[" << u1 << "]" << " <-> " << n2 << "[" << u2 << "]" << " : " << f_n << std::endl; ierr = VecSetValueLocal(f_int, n1, -f_n, ADD_VALUES); ierr = VecSetValueLocal(f_int, n2, f_n, ADD_VALUES); } ierr = VecAssemblyBegin(f_int); ierr = VecAssemblyEnd(f_int); // this->getDOFManager().assembleElementalArrayLocalArray( // forces_internal_el, internal_forces, _segment_2, ghost_type); } Real getPotentialEnergy() { std::cout << "getPotentialEnergy" << std::endl; copy(x, displacement); Vec Ax; ierr = VecDuplicate(x, &Ax); ierr = MatMult(K, x, Ax); PetscScalar res; ierr = VecDot(x, Ax, &res); return res / 2.; } Real getKineticEnergy() { std::cout << "getKineticEnergy" << std::endl; return 0; } // Real getExternalWorkIncrement() { // Real res = 0; // auto it = velocity.begin(); // auto end = velocity.end(); // auto if_it = internal_forces.begin(); // auto ef_it = forces.begin(); // auto b_it = blocked.begin(); // for (UInt node = 0; it != end; ++it, ++if_it, ++ef_it, ++b_it, ++node) { // if (mesh.isLocalOrMasterNode(node)) // res += (*b_it ? -*if_it : *ef_it) * *it; // } // mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); // return res * this->getTimeStep(); // } // void predictor() {} // void corrector() {} // /* ------------------------------------------------------------------------ // */ UInt getNbData(const Array & elements, // const SynchronizationTag &) const { // return elements.size() * sizeof(Real); // } // void packData(CommunicationBuffer & buffer, const Array & // elements, // const SynchronizationTag & tag) const { // if (tag == SynchronizationTag::_user_1) { // for (const auto & el : elements) { // buffer << this->stresses(el.element); // } // } // } // void unpackData(CommunicationBuffer & buffer, const Array & // elements, // const SynchronizationTag & tag) { // if (tag == SynchronizationTag::_user_1) { // auto cit = this->mesh.getConnectivity(_segment_2, _ghost).begin(2); // for (const auto & el : elements) { // Real stress; // buffer >> stress; // Real f = A * stress; // Vector conn = cit[el.element]; // this->internal_forces(conn(0), _x) += -f; // this->internal_forces(conn(1), _x) += f; // } // } // } Real getExternalWorkIncrement() { std::cout << "getExternalWorkIncrement" << std::endl; return 0.; } template void applyBC(Functor && func, const ID & group_id) { auto & group = mesh.getElementGroup(group_id).getNodeGroup().getNodes(); auto blocked_dofs = make_view(blocked, 1).begin(); auto disps = make_view(displacement, 1).begin(); auto poss = make_view(mesh.getNodes(), 1).begin(); for (auto && node : group) { auto disp = Vector(disps[node]); auto pos = Vector(poss[node]); auto flags = Vector(blocked_dofs[node]); func(node, flags, disp, pos); } } const Mesh & getMesh() const { return mesh; } UInt getSpatialDimension() const { return 1; } auto & getBlockedDOFs() { return blocked; } void setTimeStep(Real dt) { std::cout << "setTimeStep" << std::endl; this->dt = dt; } private: PetscErrorCode ierr{0}; MPI_Comm mpi_comm; ISLocalToGlobalMapping petsc_local_to_global; UInt nb_dofs; UInt nb_elements; bool lumped; bool is_stiffness_assembled{false}; bool is_mass_assembled{false}; bool is_lumped_mass_assembled{false}; Mat K{nullptr}, J{nullptr}, M{nullptr}; Vec rhs{nullptr}, x{nullptr}, x_save{nullptr}, dx{nullptr}, f_int{nullptr}, f_dirichlet{nullptr}; SNES snes; Real dt{0}; Array save_displacement; public: Real E, A, rho; Mesh & mesh; Array displacement; Array velocity; Array acceleration; Array blocked; Array forces; Array internal_forces; Array stresses; Array strains; Array initial_lengths; }; /* -------------------------------------------------------------------------- */ class Sinusoidal : public BC::Dirichlet::DirichletFunctor { public: Sinusoidal(MyModel & model, Real amplitude, Real pulse_width, Real t) : model(model), A(amplitude), k(2 * M_PI / pulse_width), t(t), v{std::sqrt(model.E / model.rho)} {} void operator()(UInt n, Vector & /*flags*/, Vector & disp, const Vector & coord) const { auto x = coord(_x); model.velocity(n, _x) = k * v * A * sin(k * (x - v * t)); disp(_x) = A * cos(k * (x - v * t)); } private: MyModel & model; Real A{1.}; Real k{2 * M_PI}; Real t{1.}; Real v{1.}; }; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); PetscInitialize(&argc, &argv, nullptr, nullptr); UInt prank = Communicator::getStaticCommunicator().whoAmI(); UInt global_nb_nodes = 3; UInt max_steps = 400; Real time_step = 0.001; Mesh mesh(1); Real F = -9.81; bool _explicit = EXPLICIT; const Real pulse_width = 0.2; const Real A = 0.01; if (prank == 0) genMesh(mesh, global_nb_nodes); mesh.distribute(); // mesh.makePeriodic(_x); MyModel model(F, mesh, _explicit); // model.forces.clear(); // model.blocked.clear(); // model.applyBC(Sinusoidal(model, A, pulse_width, 0.), "all"); // model.applyBC(BC::Dirichlet::FlagOnly(_x), "border"); // if (!_explicit) { // model.getNewSolver("dynamic", TimeStepSolverType::_dynamic, // NonLinearSolverType::_newton_raphson); // model.setIntegrationScheme("dynamic", "disp", // IntegrationSchemeType::_trapezoidal_rule_2, // IntegrationScheme::_displacement); // } else { // model.getNewSolver("dynamic", TimeStepSolverType::_dynamic_lumped, // NonLinearSolverType::_lumped); // model.setIntegrationScheme("dynamic", "disp", // IntegrationSchemeType::_central_difference, // IntegrationScheme::_acceleration); // } model.setTimeStep(time_step); if (prank == 0) { std::cout << std::scientific; std::cout << std::setw(14) << "time" << "," << std::setw(14) << "wext" << "," << std::setw(14) << "epot" << "," << std::setw(14) << "ekin" << "," << std::setw(14) << "total" << "," << std::setw(14) << "max_disp" << "," << std::setw(14) << "min_disp" << std::endl; } Real wext = 0.; // model.getDOFManager().clearResidual(); // model.assembleResidual(); Real epot = 0; // model.getPotentialEnergy(); Real ekin = 0; // model.getKineticEnergy(); Real einit = ekin + epot; Real etot = ekin + epot - wext - einit; Real max_disp = 0., min_disp = 0.; for (auto && disp : model.displacement) { max_disp = std::max(max_disp, disp); min_disp = std::min(min_disp, disp); } if (prank == 0) { std::cout << std::setw(14) << 0. << "," << std::setw(14) << wext << "," << std::setw(14) << epot << "," << std::setw(14) << ekin << "," << std::setw(14) << etot << "," << std::setw(14) << max_disp << "," << std::setw(14) << min_disp << std::endl; } // #if EXPLICIT == false // NonLinearSolver & solver = // model.getDOFManager().getNonLinearSolver("dynamic"); // solver.set("max_iterations", 20); // #endif auto * dumper = new DumperParaview("dynamic", "./paraview"); mesh.registerExternalDumper(*dumper, "dynamic", true); mesh.addDumpMesh(mesh); mesh.addDumpFieldExternalToDumper("dynamic", "displacement", model.displacement); mesh.addDumpFieldExternalToDumper("dynamic", "velocity", model.velocity); mesh.addDumpFieldExternalToDumper("dynamic", "forces", model.forces); mesh.addDumpFieldExternalToDumper("dynamic", "acceleration", model.acceleration); mesh.dump(); max_steps = 1; for (UInt i = 1; i < max_steps + 1; ++i) { // model.applyBC(Sinusoidal(model, A, pulse_width, time_step * (i - 1)), // "border"); model.solveStep(); mesh.dump(); epot = model.getPotentialEnergy(); ekin = model.getKineticEnergy(); wext += model.getExternalWorkIncrement(); etot = ekin + epot - wext - einit; Real max_disp = 0., min_disp = 0.; for (auto && disp : model.displacement) { max_disp = std::max(max_disp, disp); min_disp = std::min(min_disp, disp); } if (prank == 0) { std::cout << std::setw(14) << time_step * i << "," << std::setw(14) << wext << "," << std::setw(14) << epot << "," << std::setw(14) << ekin << "," << std::setw(14) << etot << "," << std::setw(14) << max_disp << "," << std::setw(14) << min_disp << std::endl; } } // output.close(); // finalize(); // PetscFinalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes) { MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); Array & conn = mesh_accessor.getConnectivity(_segment_2); nodes.resize(nb_nodes); - auto & all = mesh.createNodeGroup("all_nodes"); + //auto & all = mesh.createNodeGroup("all_nodes"); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); - all.add(n); + //all.add(n); } - mesh.createElementGroupFromNodeGroup("all", "all_nodes"); + //mesh.createElementGroupFromNodeGroup("all", "all_nodes"); conn.resize(nb_nodes - 1); for (UInt n = 0; n < nb_nodes - 1; ++n) { conn(n, 0) = n; conn(n, 1) = n + 1; } - Array & conn_points = mesh_accessor.getConnectivity(_point_1); - conn_points.resize(2); + // Array & conn_points = mesh_accessor.getConnectivity(_point_1); + // conn_points.resize(2); - conn_points(0, 0) = 0; - conn_points(1, 0) = nb_nodes - 1; + // conn_points(0, 0) = 0; + // conn_points(1, 0) = nb_nodes - 1; - auto & border = mesh.createElementGroup("border", 0); - border.add({_point_1, 0, _not_ghost}, true); - border.add({_point_1, 1, _not_ghost}, true); + // auto & border = mesh.createElementGroup("border", 0); + // border.add({_point_1, 0, _not_ghost}, true); + // border.add({_point_1, 1, _not_ghost}, true); mesh_accessor.makeReady(); } diff --git a/test/test_model/test_model_solver/test_model_solver_my_model.hh b/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh similarity index 100% rename from test/test_model/test_model_solver/test_model_solver_my_model.hh rename to test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh diff --git a/test/test_model/test_non_local_toolbox/CMakeLists.txt b/test/test_model/test_common/test_non_local_toolbox/CMakeLists.txt similarity index 100% rename from test/test_model/test_non_local_toolbox/CMakeLists.txt rename to test/test_model/test_common/test_non_local_toolbox/CMakeLists.txt diff --git a/test/test_model/test_non_local_toolbox/fine_mesh.msh b/test/test_model/test_common/test_non_local_toolbox/fine_mesh.msh similarity index 100% rename from test/test_model/test_non_local_toolbox/fine_mesh.msh rename to test/test_model/test_common/test_non_local_toolbox/fine_mesh.msh diff --git a/test/test_model/test_non_local_toolbox/material.dat b/test/test_model/test_common/test_non_local_toolbox/material.dat similarity index 100% rename from test/test_model/test_non_local_toolbox/material.dat rename to test/test_model/test_common/test_non_local_toolbox/material.dat diff --git a/test/test_model/test_non_local_toolbox/material_avg.dat b/test/test_model/test_common/test_non_local_toolbox/material_avg.dat similarity index 100% rename from test/test_model/test_non_local_toolbox/material_avg.dat rename to test/test_model/test_common/test_non_local_toolbox/material_avg.dat diff --git a/test/test_model/test_non_local_toolbox/material_parallel_test.dat b/test/test_model/test_common/test_non_local_toolbox/material_parallel_test.dat similarity index 100% rename from test/test_model/test_non_local_toolbox/material_parallel_test.dat rename to test/test_model/test_common/test_non_local_toolbox/material_parallel_test.dat diff --git a/test/test_model/test_non_local_toolbox/material_remove_damage.dat b/test/test_model/test_common/test_non_local_toolbox/material_remove_damage.dat similarity index 100% rename from test/test_model/test_non_local_toolbox/material_remove_damage.dat rename to test/test_model/test_common/test_non_local_toolbox/material_remove_damage.dat diff --git a/test/test_model/test_non_local_toolbox/material_weight_computation.dat b/test/test_model/test_common/test_non_local_toolbox/material_weight_computation.dat similarity index 100% rename from test/test_model/test_non_local_toolbox/material_weight_computation.dat rename to test/test_model/test_common/test_non_local_toolbox/material_weight_computation.dat diff --git a/test/test_model/test_non_local_toolbox/my_model.hh b/test/test_model/test_common/test_non_local_toolbox/my_model.hh similarity index 100% rename from test/test_model/test_non_local_toolbox/my_model.hh rename to test/test_model/test_common/test_non_local_toolbox/my_model.hh diff --git a/test/test_model/test_non_local_toolbox/pair_test.msh b/test/test_model/test_common/test_non_local_toolbox/pair_test.msh similarity index 100% rename from test/test_model/test_non_local_toolbox/pair_test.msh rename to test/test_model/test_common/test_non_local_toolbox/pair_test.msh diff --git a/test/test_model/test_non_local_toolbox/parallel_test.msh b/test/test_model/test_common/test_non_local_toolbox/parallel_test.msh similarity index 100% rename from test/test_model/test_non_local_toolbox/parallel_test.msh rename to test/test_model/test_common/test_non_local_toolbox/parallel_test.msh diff --git a/test/test_model/test_non_local_toolbox/plate.geo b/test/test_model/test_common/test_non_local_toolbox/plate.geo similarity index 100% rename from test/test_model/test_non_local_toolbox/plate.geo rename to test/test_model/test_common/test_non_local_toolbox/plate.geo diff --git a/test/test_model/test_non_local_toolbox/plate.msh b/test/test_model/test_common/test_non_local_toolbox/plate.msh similarity index 100% rename from test/test_model/test_non_local_toolbox/plate.msh rename to test/test_model/test_common/test_non_local_toolbox/plate.msh diff --git a/test/test_model/test_non_local_toolbox/plot_neighborhoods.py b/test/test_model/test_common/test_non_local_toolbox/plot_neighborhoods.py similarity index 100% rename from test/test_model/test_non_local_toolbox/plot_neighborhoods.py rename to test/test_model/test_common/test_non_local_toolbox/plot_neighborhoods.py diff --git a/test/test_model/test_non_local_toolbox/test_build_neighborhood_parallel.cc b/test/test_model/test_common/test_non_local_toolbox/test_build_neighborhood_parallel.cc similarity index 100% rename from test/test_model/test_non_local_toolbox/test_build_neighborhood_parallel.cc rename to test/test_model/test_common/test_non_local_toolbox/test_build_neighborhood_parallel.cc diff --git a/test/test_model/test_non_local_toolbox/test_material.cc b/test/test_model/test_common/test_non_local_toolbox/test_material.cc similarity index 100% rename from test/test_model/test_non_local_toolbox/test_material.cc rename to test/test_model/test_common/test_non_local_toolbox/test_material.cc diff --git a/test/test_model/test_non_local_toolbox/test_material.hh b/test/test_model/test_common/test_non_local_toolbox/test_material.hh similarity index 100% rename from test/test_model/test_non_local_toolbox/test_material.hh rename to test/test_model/test_common/test_non_local_toolbox/test_material.hh diff --git a/test/test_model/test_non_local_toolbox/test_material_damage.cc b/test/test_model/test_common/test_non_local_toolbox/test_material_damage.cc similarity index 100% rename from test/test_model/test_non_local_toolbox/test_material_damage.cc rename to test/test_model/test_common/test_non_local_toolbox/test_material_damage.cc diff --git a/test/test_model/test_non_local_toolbox/test_material_damage.hh b/test/test_model/test_common/test_non_local_toolbox/test_material_damage.hh similarity index 100% rename from test/test_model/test_non_local_toolbox/test_material_damage.hh rename to test/test_model/test_common/test_non_local_toolbox/test_material_damage.hh diff --git a/test/test_model/test_non_local_toolbox/test_non_local_averaging.cc b/test/test_model/test_common/test_non_local_toolbox/test_non_local_averaging.cc similarity index 100% rename from test/test_model/test_non_local_toolbox/test_non_local_averaging.cc rename to test/test_model/test_common/test_non_local_toolbox/test_non_local_averaging.cc diff --git a/test/test_model/test_non_local_toolbox/test_non_local_neighborhood_base.cc b/test/test_model/test_common/test_non_local_toolbox/test_non_local_neighborhood_base.cc similarity index 100% rename from test/test_model/test_non_local_toolbox/test_non_local_neighborhood_base.cc rename to test/test_model/test_common/test_non_local_toolbox/test_non_local_neighborhood_base.cc diff --git a/test/test_model/test_non_local_toolbox/test_non_local_neighborhood_base.verified b/test/test_model/test_common/test_non_local_toolbox/test_non_local_neighborhood_base.verified similarity index 100% rename from test/test_model/test_non_local_toolbox/test_non_local_neighborhood_base.verified rename to test/test_model/test_common/test_non_local_toolbox/test_non_local_neighborhood_base.verified diff --git a/test/test_model/test_non_local_toolbox/test_pair_computation.cc b/test/test_model/test_common/test_non_local_toolbox/test_pair_computation.cc similarity index 100% rename from test/test_model/test_non_local_toolbox/test_pair_computation.cc rename to test/test_model/test_common/test_non_local_toolbox/test_pair_computation.cc diff --git a/test/test_model/test_non_local_toolbox/test_remove_damage_weight_function.cc b/test/test_model/test_common/test_non_local_toolbox/test_remove_damage_weight_function.cc similarity index 100% rename from test/test_model/test_non_local_toolbox/test_remove_damage_weight_function.cc rename to test/test_model/test_common/test_non_local_toolbox/test_remove_damage_weight_function.cc diff --git a/test/test_model/test_non_local_toolbox/test_remove_damage_weight_function.verified b/test/test_model/test_common/test_non_local_toolbox/test_remove_damage_weight_function.verified similarity index 100% rename from test/test_model/test_non_local_toolbox/test_remove_damage_weight_function.verified rename to test/test_model/test_common/test_non_local_toolbox/test_remove_damage_weight_function.verified diff --git a/test/test_model/test_non_local_toolbox/test_weight_computation.cc b/test/test_model/test_common/test_non_local_toolbox/test_weight_computation.cc similarity index 100% rename from test/test_model/test_non_local_toolbox/test_weight_computation.cc rename to test/test_model/test_common/test_non_local_toolbox/test_weight_computation.cc diff --git a/test/test_model/test_non_local_toolbox/test_weight_computation.verified b/test/test_model/test_common/test_non_local_toolbox/test_weight_computation.verified similarity index 100% rename from test/test_model/test_non_local_toolbox/test_weight_computation.verified rename to test/test_model/test_common/test_non_local_toolbox/test_weight_computation.verified