diff --git a/src/model/common/dof_manager/dof_manager_default.cc b/src/model/common/dof_manager/dof_manager_default.cc index 08412cb48..c74669ba4 100644 --- a/src/model/common/dof_manager/dof_manager_default.cc +++ b/src/model/common/dof_manager/dof_manager_default.cc @@ -1,476 +1,485 @@ /** * @file dof_manager_default.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Thu Feb 08 2018 * * @brief Implementation of the default DOFManager * * @section LICENSE * * Copyright (©) 2015-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 "dof_manager_default.hh" #include "communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include "non_linear_solver_default.hh" #include "periodic_node_synchronizer.hh" #include "solver_vector_default.hh" #include "solver_vector_distributed.hh" #include "sparse_matrix_aij.hh" #include "time_step_solver_default.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(const ID & id, const MemoryID & memory_id) : DOFManager(id, memory_id), synchronizer(nullptr) { residual = std::make_unique( *this, std::string(id + ":residual")); solution = std::make_unique( *this, std::string(id + ":solution")); data_cache = std::make_unique( *this, std::string(id + ":data_cache")); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(Mesh & mesh, const ID & id, const MemoryID & memory_id) : DOFManager(mesh, id, memory_id), synchronizer(nullptr) { if (this->mesh->isDistributed()) { this->synchronizer = std::make_unique( *this, this->id + ":dof_synchronizer", this->memory_id); residual = std::make_unique( *this, std::string(id + ":residual")); solution = std::make_unique( *this, std::string(id + ":solution")); data_cache = std::make_unique( *this, std::string(id + ":data_cache")); } else { residual = std::make_unique( *this, std::string(id + ":residual")); solution = std::make_unique( *this, std::string(id + ":solution")); data_cache = std::make_unique( *this, std::string(id + ":data_cache")); } } /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() = default; /* -------------------------------------------------------------------------- */ void DOFManagerDefault::makeConsistentForPeriodicity(const ID & dof_id, SolverVector & array) { auto & dof_data = this->getDOFDataTyped(dof_id); if (dof_data.support_type != _dst_nodal) return; if (not mesh->isPeriodic()) return; this->mesh->getPeriodicNodeSynchronizer() .reduceSynchronizeWithPBCSlaves( aka::as_type(array).getVector()); } /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor) { AKANTU_DEBUG_IN(); auto & dof_data = this->getDOFDataTyped(dof_id); AKANTU_DEBUG_ASSERT(dof_data.local_equation_number.size() == array_to_assemble.size() * array_to_assemble.getNbComponent(), "The array to assemble does not have a correct size." << " (" << array_to_assemble.getID() << ")"); if (dof_data.support_type == _dst_nodal and mesh->isPeriodic()) { for (auto && data : zip(dof_data.local_equation_number, dof_data.associated_nodes, make_view(array_to_assemble))) { auto && equ_num = std::get<0>(data); auto && node = std::get<1>(data); auto && arr = std::get<2>(data); global_array(equ_num) += scale_factor * (arr) * (not this->mesh->isPeriodicSlave(node)); } } else { for (auto && data : zip(dof_data.local_equation_number, make_view(array_to_assemble))) { auto && equ_num = std::get<0>(data); auto && arr = std::get<1>(data); global_array(equ_num) += scale_factor * (arr); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, SolverVector & global_array_v, Real scale_factor) { assembleToGlobalArray( dof_id, array_to_assemble, aka::as_type(global_array_v).getVector(), scale_factor); } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFDataDefault::DOFDataDefault(const ID & dof_id) : DOFData(dof_id) {} /* -------------------------------------------------------------------------- */ auto DOFManagerDefault::getNewDOFData(const ID & dof_id) -> std::unique_ptr { return std::make_unique(dof_id); } /* -------------------------------------------------------------------------- */ std::tuple DOFManagerDefault::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { auto ret = DOFManager::registerDOFsInternal(dof_id, dofs_array); // update the synchronizer if needed if (this->synchronizer) this->synchronizer->registerDOFs(dof_id); return ret; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const MatrixType & matrix_type) { return this->registerSparseMatrix(*this, id, matrix_type); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const ID & matrix_to_copy_id) { return this->registerSparseMatrix(id, matrix_to_copy_id); } /* -------------------------------------------------------------------------- */ SolverVector & DOFManagerDefault::getNewLumpedMatrix(const ID & id) { return this->registerLumpedMatrix(*this, id); } /* -------------------------------------------------------------------------- */ SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & id) { auto & matrix = DOFManager::getMatrix(id); return aka::as_type(matrix); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManagerDefault::getNewNonLinearSolver(const ID & id, const NonLinearSolverType & type) { switch (type) { #if defined(AKANTU_IMPLICIT) case NonLinearSolverType::_newton_raphson: /* FALLTHRU */ /* [[fallthrough]]; un-comment when compiler will get it */ case NonLinearSolverType::_newton_raphson_modified: { return this->registerNonLinearSolver( *this, id, type); } case NonLinearSolverType::_linear: { return this->registerNonLinearSolver(*this, id, type); } #endif case NonLinearSolverType::_lumped: { return this->registerNonLinearSolver(*this, id, type); } default: AKANTU_EXCEPTION("The asked type of non linear solver is not supported by " "this dof manager"); } } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { return this->registerTimeStepSolver(*this, id, type, non_linear_solver); } /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id, const Array & global_array, Array & local_array) const { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationsNumbers(dof_id); UInt nb_degree_of_freedoms = equation_number.size(); local_array.resize(nb_degree_of_freedoms / local_array.getNbComponent()); auto loc_it = local_array.begin_reinterpret(nb_degree_of_freedoms); auto equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++loc_it, ++equ_it) { (*loc_it) = global_array(*equ_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getArrayPerDOFs(const ID & dof_id, const SolverVector & global_array, Array & local_array) { getArrayPerDOFs( dof_id, aka::as_type(global_array).getVector(), local_array); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleLumpedMatMulVectToResidual( const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { const Array & A = this->getLumpedMatrix(A_id); auto & cache = aka::as_type(*this->data_cache); cache.clear(); this->assembleToGlobalArray(dof_id, x, cache.getVector(), scale_factor); for (auto && data : zip(make_view(A), make_view(cache.getVector()), make_view(this->getResidualArray()))) { const auto & A = std::get<0>(data); const auto & x = std::get<1>(data); auto & r = std::get<2>(data); r += A * x; } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { this->addToProfile(matrix_id, dof_id, type, ghost_type); auto & A = getMatrix(matrix_id); DOFManager::assembleElementalMatricesToMatrix_( A, dof_id, elementary_mat, type, ghost_type, elemental_matrix_type, filter_elements); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assemblePreassembledMatrix( const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) { auto & A = getMatrix(matrix_id); DOFManager::assemblePreassembledMatrix_(A, dof_id_m, dof_id_n, terms); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleMatMulVectToArray(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor) { if (mesh->isDistributed()) { DOFManager::assembleMatMulVectToArray_( dof_id, A_id, x, array, scale_factor); } else { DOFManager::assembleMatMulVectToArray_( dof_id, A_id, x, array, scale_factor); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const auto & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) return; auto mat_dof = std::make_pair(matrix_id, dof_id); auto type_pair = std::make_pair(type, ghost_type); auto prof_it = this->matrix_profiled_dofs.find(mat_dof); if (prof_it != this->matrix_profiled_dofs.end() && std::find(prof_it->second.begin(), prof_it->second.end(), type_pair) != prof_it->second.end()) return; auto nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent(); const auto & equation_number = this->getLocalEquationsNumbers(dof_id); auto & A = this->getMatrix(matrix_id); auto size = A.size(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const auto & connectivity = this->mesh->getConnectivity(type, ghost_type); auto cbegin = connectivity.begin(nb_nodes_per_element); auto cit = cbegin; auto nb_elements = connectivity.size(); UInt * ge_it = nullptr; if (dof_data.group_support != "__mesh__") { const auto & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); ge_it = group_elements.storage(); nb_elements = group_elements.size(); } UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node; Vector element_eq_nb(size_mat); for (UInt e = 0; e < nb_elements; ++e) { if (ge_it) cit = cbegin + *ge_it; this->extractElementEquationNumber( equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb); std::transform( element_eq_nb.storage(), element_eq_nb.storage() + element_eq_nb.size(), element_eq_nb.storage(), [&](auto & local) { return this->localToGlobalEquationNumber(local); }); if (ge_it) ++ge_it; else ++cit; for (UInt i = 0; i < size_mat; ++i) { UInt c_irn = element_eq_nb(i); if (c_irn < size) { for (UInt j = 0; j < size_mat; ++j) { UInt c_jcn = element_eq_nb(j); if (c_jcn < size) { A.add(c_irn, c_jcn); } } } } } this->matrix_profiled_dofs[mat_dof].push_back(type_pair); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Array & DOFManagerDefault::getSolutionArray() { return dynamic_cast(this->solution.get())->getVector(); } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getResidualArray() const { return dynamic_cast(this->residual.get())->getVector(); } /* -------------------------------------------------------------------------- */ Array & DOFManagerDefault::getResidualArray() { return dynamic_cast(this->residual.get())->getVector(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { DOFManager::onNodesAdded(nodes_list, event); if (this->synchronizer) this->synchronizer->onNodesAdded(nodes_list); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::resizeGlobalArrays() { DOFManager::resizeGlobalArrays(); this->global_blocked_dofs.resize(this->local_system_size, true); this->previous_global_blocked_dofs.resize(this->local_system_size, true); matrix_profiled_dofs.clear(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateGlobalBlockedDofs() { DOFManager::updateGlobalBlockedDofs(); if (this->global_blocked_dofs_release == this->previous_global_blocked_dofs_release) return; global_blocked_dofs_uint.resize(local_system_size); global_blocked_dofs_uint.set(false); for (const auto & dof : global_blocked_dofs) { global_blocked_dofs_uint[dof] = true; } } /* -------------------------------------------------------------------------- */ Array & DOFManagerDefault::getBlockedDOFs() { return global_blocked_dofs_uint; } /* -------------------------------------------------------------------------- */ const Array & DOFManagerDefault::getBlockedDOFs() const { return global_blocked_dofs_uint; } /* -------------------------------------------------------------------------- */ // register in factory -static bool default_dof_manager_is_registered [[gnu::unused]] = - DefaultDOFManagerFactory::getInstance().registerAllocator( +// static bool default_dof_manager_is_registered [[gnu::unused]] = +// DefaultDOFManagerFactory::getInstance().registerAllocator( +// "default", +// [](const ID & id, +// const MemoryID & mem_id) -> std::unique_ptr { +// return std::make_unique(id, mem_id); +// }); + +static bool dof_manager_is_registered [[gnu::unused]] = + DOFManagerFactory::getInstance().registerAllocator( "default", - [](const ID & id, + [](Mesh & mesh, const ID & id, const MemoryID & mem_id) -> std::unique_ptr { - return std::make_unique(id, mem_id); + return std::make_unique(mesh, id, mem_id); }); -static bool dof_manager_is_registered [[gnu::unused]] = +static bool dof_manager_is_registered_mumps [[gnu::unused]] = DOFManagerFactory::getInstance().registerAllocator( - "default", + "mumps", [](Mesh & mesh, const ID & id, const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(mesh, id, mem_id); }); + } // namespace akantu diff --git a/src/model/common/dof_manager/dof_manager_petsc.cc b/src/model/common/dof_manager/dof_manager_petsc.cc index 47a1f4986..f5dba87ae 100644 --- a/src/model/common/dof_manager/dof_manager_petsc.cc +++ b/src/model/common/dof_manager/dof_manager_petsc.cc @@ -1,316 +1,313 @@ /** * @file dof_manager_petsc.cc * * @author Nicolas Richart * * @date creation: Wed Oct 07 2015 * @date last modification: Tue Feb 20 2018 * * @brief DOFManaterPETSc is the PETSc implementation of the DOFManager * * @section LICENSE * * Copyright (©) 2015-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 "dof_manager_petsc.hh" #include "aka_iterators.hh" #include "communicator.hh" #include "cppargparse.hh" #include "non_linear_solver_petsc.hh" #include "solver_vector_petsc.hh" #include "sparse_matrix_petsc.hh" #include "time_step_solver_default.hh" #if defined(AKANTU_USE_MPI) #include "mpi_communicator_data.hh" #endif /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { class PETScSingleton { private: PETScSingleton() { PETSc_call(PetscInitialized, &is_initialized); if (not is_initialized) { cppargparse::ArgumentParser & argparser = getStaticArgumentParser(); int & argc = argparser.getArgC(); char **& argv = argparser.getArgV(); PETSc_call(PetscInitialize, &argc, &argv, NULL, NULL); PETSc_call( PetscPopErrorHandler); // remove the default PETSc signal handler PETSc_call(PetscPushErrorHandler, PetscIgnoreErrorHandler, NULL); } } public: PETScSingleton(const PETScSingleton &) = delete; PETScSingleton & operator=(const PETScSingleton &) = delete; ~PETScSingleton() { if (not is_initialized) { PetscFinalize(); } } static PETScSingleton & getInstance() { static PETScSingleton instance; return instance; } private: PetscBool is_initialized; }; /* -------------------------------------------------------------------------- */ DOFManagerPETSc::DOFDataPETSc::DOFDataPETSc(const ID & dof_id) : DOFData(dof_id) {} /* -------------------------------------------------------------------------- */ DOFManagerPETSc::DOFManagerPETSc(const ID & id, const MemoryID & memory_id) : DOFManager(id, memory_id) { init(); } /* -------------------------------------------------------------------------- */ DOFManagerPETSc::DOFManagerPETSc(Mesh & mesh, const ID & id, const MemoryID & memory_id) : DOFManager(mesh, id, memory_id) { init(); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::init() { // check if the akantu types and PETSc one are consistant static_assert(sizeof(Int) == sizeof(PetscInt), "The integer type of Akantu does not match the one from PETSc"); static_assert(sizeof(Real) == sizeof(PetscReal), "The integer type of Akantu does not match the one from PETSc"); #if defined(AKANTU_USE_MPI) const auto & mpi_data = aka::as_type(communicator.getCommunicatorData()); MPI_Comm mpi_comm = mpi_data.getMPICommunicator(); this->mpi_communicator = mpi_comm; #else this->mpi_communicator = PETSC_COMM_SELF; #endif PETScSingleton & instance [[gnu::unused]] = PETScSingleton::getInstance(); } /* -------------------------------------------------------------------------- */ DOFManagerPETSc::~DOFManagerPETSc() { // if (is_ltog_map) // PETSc_call(ISLocalToGlobalMappingDestroy, &is_ltog_map); } /* -------------------------------------------------------------------------- */ auto DOFManagerPETSc::getNewDOFData(const ID & dof_id) -> std::unique_ptr { return std::make_unique(dof_id); } /* -------------------------------------------------------------------------- */ std::tuple DOFManagerPETSc::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { dofs_ids.push_back(dof_id); auto ret = DOFManager::registerDOFsInternal(dof_id, dofs_array); UInt nb_dofs, nb_pure_local_dofs; std::tie(nb_dofs, nb_pure_local_dofs, std::ignore) = ret; auto && vector = std::make_unique(*this, id + ":solution"); auto x = vector->getVec(); PETSc_call(VecGetLocalToGlobalMapping, x, &is_ltog_map); // redoing the indexes based on the petsc numbering for (auto & dof_id : dofs_ids) { auto & dof_data = this->getDOFDataTyped(dof_id); Array gidx(dof_data.local_equation_number.size()); for (auto && data : zip(dof_data.local_equation_number, gidx)) { std::get<1>(data) = localToGlobalEquationNumber(std::get<0>(data)); } auto & lidx = dof_data.local_equation_number_petsc; if (is_ltog_map) { lidx.resize(gidx.size()); PetscInt n; PETSc_call(ISGlobalToLocalMappingApply, is_ltog_map, IS_GTOLM_MASK, gidx.size(), gidx.storage(), &n, lidx.storage()); } } residual = std::make_unique(*vector, id + ":residual"); data_cache = std::make_unique(*vector, id + ":data_cache"); solution = std::move(vector); for(auto & mat : matrices) { auto & A = this->getMatrix(mat.first); A.resize(); } - std::cout << "Allocating vectors " << solution.get() << std::endl; - // should also redo the lumped matrix - return ret; } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, SolverVector & global_array, Real scale_factor) { const auto & dof_data = getDOFDataTyped(dof_id); auto & g = aka::as_type(global_array); AKANTU_DEBUG_ASSERT(array_to_assemble.size() * array_to_assemble.getNbComponent() == dof_data.local_nb_dofs, "The array to assemble does not have the proper size"); g.addValuesLocal(dof_data.local_equation_number_petsc, array_to_assemble, scale_factor); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::getArrayPerDOFs(const ID & dof_id, const SolverVector & global_array, Array & local) { const auto & dof_data = getDOFDataTyped(dof_id); const auto & petsc_vector = aka::as_type(global_array); AKANTU_DEBUG_ASSERT( local.size() * local.getNbComponent() == dof_data.local_nb_dofs, "The array to get the values does not have the proper size"); petsc_vector.getValuesLocal(dof_data.local_equation_number_petsc, local); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { auto & A = getMatrix(matrix_id); DOFManager::assembleElementalMatricesToMatrix_( A, dof_id, elementary_mat, type, ghost_type, elemental_matrix_type, filter_elements); A.applyModifications(); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::assemblePreassembledMatrix( const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) { auto & A = getMatrix(matrix_id); DOFManager::assemblePreassembledMatrix_(A, dof_id_m, dof_id_n, terms); A.applyModifications(); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::assembleMatMulVectToArray(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor) { DOFManager::assembleMatMulVectToArray_( dof_id, A_id, x, array, scale_factor); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::makeConsistentForPeriodicity(const ID & /*dof_id*/, SolverVector & /*array*/) {} /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManagerPETSc::getNewNonLinearSolver(const ID & id, const NonLinearSolverType & type) { return this->registerNonLinearSolver(*this, id, type); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerPETSc::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { return this->registerTimeStepSolver(*this, id, type, non_linear_solver); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerPETSc::getNewMatrix(const ID & id, const MatrixType & matrix_type) { return this->registerSparseMatrix(*this, id, matrix_type); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerPETSc::getNewMatrix(const ID & id, const ID & matrix_to_copy_id) { return this->registerSparseMatrix(id, matrix_to_copy_id); } /* -------------------------------------------------------------------------- */ SparseMatrixPETSc & DOFManagerPETSc::getMatrix(const ID & id) { auto & matrix = DOFManager::getMatrix(id); return aka::as_type(matrix); } /* -------------------------------------------------------------------------- */ SolverVector & DOFManagerPETSc::getNewLumpedMatrix(const ID & id) { return this->registerLumpedMatrix(*this, id); } /* -------------------------------------------------------------------------- */ SolverVectorPETSc & DOFManagerPETSc::getSolution() { return aka::as_type(*this->solution); } const SolverVectorPETSc & DOFManagerPETSc::getSolution() const { return aka::as_type(*this->solution); } SolverVectorPETSc & DOFManagerPETSc::getResidual() { return aka::as_type(*this->residual); } const SolverVectorPETSc & DOFManagerPETSc::getResidual() const { return aka::as_type(*this->residual); } /* -------------------------------------------------------------------------- */ static bool dof_manager_is_registered [[gnu::unused]] = DOFManagerFactory::getInstance().registerAllocator( "petsc", [](Mesh & mesh, const ID & id, const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(mesh, id, mem_id); }); } // namespace akantu diff --git a/test/test_model/test_common/test_model_solver/CMakeLists.txt b/test/test_model/test_common/test_model_solver/CMakeLists.txt index 0bce80cc4..52eff679c 100644 --- a/test/test_model/test_common/test_model_solver/CMakeLists.txt +++ b/test/test_model/test_common/test_model_solver/CMakeLists.txt @@ -1,55 +1,62 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date creation: Fri Sep 03 2010 # @date last modification: Sat Apr 01 2017 # # @brief test for the common solvers interface of the models # # @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 . # #=============================================================================== register_test(test_dof_manager_default SOURCES test_dof_manager_default.cc PACKAGE mumps ) -register_test(test_model_solver +register_test(test_model_solver_mumps SOURCES test_model_solver.cc + COMPILE_OPTIONS "DOF_MANAGER_TYPE=\"mumps\"" PACKAGE mumps ) +register_test(test_model_solver_petsc + SOURCES test_model_solver.cc + COMPILE_OPTIONS "DOF_MANAGER_TYPE=\"petsc\"" + PACKAGE petsc + ) + register_test(test_model_solver_dynamic_explicit SOURCES test_model_solver_dynamic.cc PACKAGE core COMPILE_OPTIONS "EXPLICIT=true" ) register_test(test_model_solver_dynamic_implicit SOURCES test_model_solver_dynamic.cc PACKAGE mumps COMPILE_OPTIONS "EXPLICIT=false" ) register_test(test_model_solver_dynamic_petsc - SOURCES test_model_solver_dynamic_petsc.cc + SOURCES test_model_solver_dynamic.cc PACKAGE petsc - COMPILE_OPTIONS "EXPLICIT=false" + COMPILE_OPTIONS "EXPLICIT=false;DOF_MANAGER_TYPE=\"petsc\"" ) register_test(test_vector_petsc SOURCES test_vector_petsc.cc PACKAGE petsc ) diff --git a/test/test_model/test_common/test_model_solver/test_model_solver.cc b/test/test_model/test_common/test_model_solver/test_model_solver.cc index ccd255b27..1db74aea2 100644 --- a/test/test_model/test_common/test_model_solver/test_model_solver.cc +++ b/test/test_model/test_common/test_model_solver/test_model_solver.cc @@ -1,170 +1,175 @@ /** * @file test_model_solver.cc * * @author Nicolas Richart * * @date creation: Wed Apr 13 2016 * @date last modification: Thu Feb 01 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 "aka_random_generator.hh" #include "dof_manager.hh" #include "dof_synchronizer.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "model_solver.hh" #include "non_linear_solver.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #include "test_model_solver_my_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; static void genMesh(Mesh & mesh, UInt nb_nodes); static void printResults(MyModel & model, UInt nb_nodes); Real F = -10; /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); UInt prank = Communicator::getStaticCommunicator().whoAmI(); std::cout << std::setprecision(7); - UInt global_nb_nodes = 5; + ID dof_manager_type = "default"; +#if defined(DOF_MANAGER_TYPE) + dof_manager_type = DOF_MANAGER_TYPE; +#endif + + UInt global_nb_nodes = 100; Mesh mesh(1); RandomGenerator::seed(1); if (prank == 0) { genMesh(mesh, global_nb_nodes); } - + // std::cout << prank << RandGenerator::seed() << std::endl; mesh.distribute(); - MyModel model(F, mesh, false); + MyModel model(F, mesh, false, dof_manager_type); model.getNewSolver("static", TimeStepSolverType::_static, NonLinearSolverType::_newton_raphson); model.setIntegrationScheme("static", "disp", IntegrationSchemeType::_pseudo_time); NonLinearSolver & solver = model.getDOFManager().getNonLinearSolver("static"); solver.set("max_iterations", 2); model.solveStep(); printResults(model, global_nb_nodes); finalize(); 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); mesh_accessor.setNbGlobalNodes(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); } conn.resize(nb_nodes - 1); for (UInt n = 0; n < nb_nodes - 1; ++n) { conn(n, 0) = n; conn(n, 1) = n + 1; } mesh_accessor.makeReady(); } /* -------------------------------------------------------------------------- */ -void printResults(MyModel & model, UInt nb_nodes) { +void printResults(MyModel & model, UInt /*nb_nodes*/) { // if (model.mesh.isDistributed()) { // UInt prank = model.mesh.getCommunicator().whoAmI(); // auto & sync = dynamic_cast(model.getDOFManager()) // .getSynchronizer(); // if (prank == 0) { // Array global_displacement(nb_nodes); // Array global_forces(nb_nodes); // Array global_blocked(nb_nodes); // sync.gather(model.forces, global_forces); // sync.gather(model.displacement, global_displacement); // sync.gather(model.blocked, global_blocked); // auto force_it = global_forces.begin(); // auto disp_it = global_displacement.begin(); // auto blocked_it = global_blocked.begin(); // std::cout << "node" // << ", " << std::setw(8) << "disp" // << ", " << std::setw(8) << "force" // << ", " << std::setw(8) << "blocked" << std::endl; // UInt node = 0; // for (; disp_it != global_displacement.end(); // ++disp_it, ++force_it, ++blocked_it, ++node) { // std::cout << node << ", " << std::setw(8) << *disp_it << ", " // << std::setw(8) << *force_it << ", " << std::setw(8) // << *blocked_it << std::endl; // std::cout << std::flush; // } // } else { // sync.gather(model.forces); // sync.gather(model.displacement); // sync.gather(model.blocked); // } // } else { auto force_it = model.forces.begin(); auto disp_it = model.displacement.begin(); auto blocked_it = model.blocked.begin(); std::cout << "node" << ", " << std::setw(8) << "disp" << ", " << std::setw(8) << "force" << ", " << std::setw(8) << "blocked" << std::endl; UInt node = 0; for (; disp_it != model.displacement.end(); ++disp_it, ++force_it, ++blocked_it, ++node) { std::cout << node << ", " << std::setw(8) << *disp_it << ", " << std::setw(8) << *force_it << ", " << std::setw(8) << *blocked_it << std::endl; std::cout << std::flush; } // } } diff --git a/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc index feb9a7319..1f41000db 100644 --- a/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc +++ b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc @@ -1,240 +1,245 @@ /** * @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 "dumpable_inline_impl.hh" #include "dumper_element_partition.hh" #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include "test_model_solver_my_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef EXPLICIT #define EXPLICIT true #endif using namespace akantu; 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.}; }; static void genMesh(Mesh & mesh, UInt nb_nodes); /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); UInt prank = Communicator::getStaticCommunicator().whoAmI(); UInt global_nb_nodes = 201; 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; + ID dof_manager_type = "default"; +#if defined(DOF_MANAGER_TYPE) + dof_manager_type = DOF_MANAGER_TYPE; +#endif + if (prank == 0) genMesh(mesh, global_nb_nodes); mesh.distribute(); // mesh.makePeriodic(_x); - MyModel model(F, mesh, _explicit); + MyModel model(F, mesh, _explicit, dof_manager_type); 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", "internal_forces", model.internal_forces); mesh.addDumpFieldExternalToDumper("dynamic", "acceleration", model.acceleration); mesh.dump(); for (UInt i = 1; i < max_steps + 1; ++i) { model.applyBC(Sinusoidal(model, A, pulse_width, time_step * (i - 1)), "border"); model.solveStep("dynamic"); 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(); 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"); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); all.add(n); } 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); 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); mesh_accessor.makeReady(); } diff --git a/test/test_model/test_common/test_model_solver/test_model_solver_dynamic_petsc.verified b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic_petsc.verified new file mode 100644 index 000000000..80183659d --- /dev/null +++ b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic_petsc.verified @@ -0,0 +1,402 @@ + time, wext, epot, ekin, total, max_disp, min_disp + 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 1.000000e-02, -1.000000e-02 + 1.000000e-03, 1.325501e-20, 2.463551e-02, 2.455451e-02, 4.919002e-02, 1.000000e-02, -9.997528e-03 + 2.000000e-03, 5.264513e-07, 2.465973e-02, 2.453049e-02, 4.918970e-02, 9.995066e-03, -9.987643e-03 + 3.000000e-03, 1.646899e-06, 2.468359e-02, 2.450696e-02, 4.918890e-02, 9.992367e-03, -9.992367e-03 + 4.000000e-03, 3.434889e-06, 2.470705e-02, 2.448383e-02, 4.918744e-02, 1.001196e-02, -1.001196e-02 + 5.000000e-03, 5.954621e-06, 2.473006e-02, 2.446103e-02, 4.918513e-02, 1.002165e-02, -1.002165e-02 + 6.000000e-03, 9.248431e-06, 2.475255e-02, 2.443853e-02, 4.918184e-02, 1.002143e-02, -1.002143e-02 + 7.000000e-03, 1.333043e-05, 2.477445e-02, 2.441638e-02, 4.917749e-02, 1.001130e-02, -1.001130e-02 + 8.000000e-03, 1.818752e-05, 2.479568e-02, 2.439463e-02, 4.917213e-02, 1.001456e-02, -1.001456e-02 + 9.000000e-03, 2.378653e-05, 2.481620e-02, 2.437340e-02, 4.916581e-02, 1.003376e-02, -1.003376e-02 + 1.000000e-02, 3.008447e-05, 2.483591e-02, 2.435282e-02, 4.915865e-02, 1.004302e-02, -1.004302e-02 + 1.100000e-02, 3.703825e-05, 2.485476e-02, 2.433301e-02, 4.915073e-02, 1.004236e-02, -1.004236e-02 + 1.200000e-02, 4.461078e-05, 2.487267e-02, 2.431407e-02, 4.914213e-02, 1.003177e-02, -1.003177e-02 + 1.300000e-02, 5.277227e-05, 2.488956e-02, 2.429609e-02, 4.913288e-02, 1.003231e-02, -1.003231e-02 + 1.400000e-02, 6.149723e-05, 2.490537e-02, 2.427910e-02, 4.912298e-02, 1.005091e-02, -1.005091e-02 + 1.500000e-02, 7.075925e-05, 2.492003e-02, 2.426315e-02, 4.911242e-02, 1.005956e-02, -1.005956e-02 + 1.600000e-02, 8.052626e-05, 2.493348e-02, 2.424827e-02, 4.910123e-02, 1.005827e-02, -1.005827e-02 + 1.700000e-02, 9.075818e-05, 2.494567e-02, 2.423452e-02, 4.908943e-02, 1.004704e-02, -1.004704e-02 + 1.800000e-02, 1.014079e-04, 2.495654e-02, 2.422198e-02, 4.907711e-02, 1.004388e-02, -1.004388e-02 + 1.900000e-02, 1.124246e-04, 2.496604e-02, 2.421072e-02, 4.906434e-02, 1.006174e-02, -1.006174e-02 + 2.000000e-02, 1.237583e-04, 2.497414e-02, 2.420082e-02, 4.905119e-02, 1.006964e-02, -1.006964e-02 + 2.100000e-02, 1.353632e-04, 2.498079e-02, 2.419233e-02, 4.903775e-02, 1.006760e-02, -1.006760e-02 + 2.200000e-02, 1.471981e-04, 2.498595e-02, 2.418529e-02, 4.902405e-02, 1.005559e-02, -1.005559e-02 + 2.300000e-02, 1.592255e-04, 2.498962e-02, 2.417972e-02, 4.901011e-02, 1.004814e-02, -1.004814e-02 + 2.400000e-02, 1.714076e-04, 2.499175e-02, 2.417563e-02, 4.899597e-02, 1.006519e-02, -1.006519e-02 + 2.500000e-02, 1.837036e-04, 2.499234e-02, 2.417300e-02, 4.898164e-02, 1.007228e-02, -1.007228e-02 + 2.600000e-02, 1.960671e-04, 2.499138e-02, 2.417186e-02, 4.896717e-02, 1.006942e-02, -1.006942e-02 + 2.700000e-02, 2.084467e-04, 2.498887e-02, 2.417221e-02, 4.895263e-02, 1.005659e-02, -1.005659e-02 + 2.800000e-02, 2.207882e-04, 2.498480e-02, 2.417408e-02, 4.893809e-02, 1.004470e-02, -1.004470e-02 + 2.900000e-02, 2.330383e-04, 2.497918e-02, 2.417748e-02, 4.892363e-02, 1.006093e-02, -1.006093e-02 + 3.000000e-02, 2.451469e-04, 2.497204e-02, 2.418243e-02, 4.890932e-02, 1.006722e-02, -1.006722e-02 + 3.100000e-02, 2.570693e-04, 2.496339e-02, 2.418891e-02, 4.889523e-02, 1.006355e-02, -1.006355e-02 + 3.200000e-02, 2.687651e-04, 2.495326e-02, 2.419689e-02, 4.888138e-02, 1.004993e-02, -1.004993e-02 + 3.300000e-02, 2.801963e-04, 2.494167e-02, 2.420633e-02, 4.886781e-02, 1.003389e-02, -1.003405e-02 + 3.400000e-02, 2.913245e-04, 2.492867e-02, 2.421719e-02, 4.885454e-02, 1.004939e-02, -1.004951e-02 + 3.500000e-02, 3.021093e-04, 2.491431e-02, 2.422941e-02, 4.884161e-02, 1.005495e-02, -1.005496e-02 + 3.600000e-02, 3.125073e-04, 2.489863e-02, 2.424295e-02, 4.882907e-02, 1.005057e-02, -1.005057e-02 + 3.700000e-02, 3.224742e-04, 2.488168e-02, 2.425778e-02, 4.881699e-02, 1.003625e-02, -1.003625e-02 + 3.800000e-02, 3.319667e-04, 2.486353e-02, 2.427386e-02, 4.880542e-02, 1.001679e-02, -1.001747e-02 + 3.900000e-02, 3.409452e-04, 2.484424e-02, 2.429114e-02, 4.879444e-02, 1.003169e-02, -1.003314e-02 + 4.000000e-02, 3.493756e-04, 2.482388e-02, 2.430957e-02, 4.878408e-02, 1.003668e-02, -1.003882e-02 + 4.100000e-02, 3.572291e-04, 2.480253e-02, 2.432907e-02, 4.877437e-02, 1.003174e-02, -1.003423e-02 + 4.200000e-02, 3.644813e-04, 2.478025e-02, 2.434957e-02, 4.876534e-02, 1.001689e-02, -1.001916e-02 + 4.300000e-02, 3.711096e-04, 2.475715e-02, 2.437097e-02, 4.875701e-02, 9.995079e-03, -9.995079e-03 + 4.400000e-02, 3.770919e-04, 2.473329e-02, 2.439318e-02, 4.874937e-02, 1.000959e-02, -1.000959e-02 + 4.500000e-02, 3.824052e-04, 2.470877e-02, 2.441611e-02, 4.874247e-02, 1.001420e-02, -1.001420e-02 + 4.600000e-02, 3.870268e-04, 2.468368e-02, 2.443969e-02, 4.873634e-02, 1.000891e-02, -1.000891e-02 + 4.700000e-02, 3.909357e-04, 2.465811e-02, 2.446382e-02, 4.873100e-02, 9.993727e-03, -9.993727e-03 + 4.800000e-02, 3.941142e-04, 2.463217e-02, 2.448845e-02, 4.872650e-02, 9.970899e-03, -9.989870e-03 + 4.900000e-02, 3.965501e-04, 2.460595e-02, 2.451347e-02, 4.872286e-02, 9.985244e-03, -1.000478e-02 + 5.000000e-02, 3.982368e-04, 2.457954e-02, 2.453879e-02, 4.872009e-02, 9.989716e-03, -1.000881e-02 + 5.100000e-02, 3.991723e-04, 2.455306e-02, 2.456430e-02, 4.871819e-02, 9.984310e-03, -1.000224e-02 + 5.200000e-02, 3.993585e-04, 2.452660e-02, 2.458990e-02, 4.871714e-02, 9.969033e-03, -9.985367e-03 + 5.300000e-02, 3.987985e-04, 2.450027e-02, 2.461547e-02, 4.871693e-02, 9.946627e-03, -9.958363e-03 + 5.400000e-02, 3.974965e-04, 2.447416e-02, 2.464091e-02, 4.871757e-02, 9.961050e-03, -9.961050e-03 + 5.500000e-02, 3.954573e-04, 2.444838e-02, 2.466612e-02, 4.871904e-02, 9.965623e-03, -9.965623e-03 + 5.600000e-02, 3.926871e-04, 2.442303e-02, 2.469101e-02, 4.872135e-02, 9.960343e-03, -9.960343e-03 + 5.700000e-02, 3.891953e-04, 2.439821e-02, 2.471548e-02, 4.872450e-02, 9.945215e-03, -9.945215e-03 + 5.800000e-02, 3.849955e-04, 2.437402e-02, 2.473946e-02, 4.872848e-02, 9.924653e-03, -9.924653e-03 + 5.900000e-02, 3.801062e-04, 2.435055e-02, 2.476284e-02, 4.873328e-02, 9.939383e-03, -9.939383e-03 + 6.000000e-02, 3.745503e-04, 2.432790e-02, 2.478552e-02, 4.873886e-02, 9.944286e-03, -9.944286e-03 + 6.100000e-02, 3.683542e-04, 2.430615e-02, 2.480741e-02, 4.874520e-02, 9.939356e-03, -9.939356e-03 + 6.200000e-02, 3.615462e-04, 2.428540e-02, 2.482840e-02, 4.875225e-02, 9.924599e-03, -9.924599e-03 + 6.300000e-02, 3.541557e-04, 2.426572e-02, 2.484842e-02, 4.875998e-02, 9.907140e-03, -9.907138e-03 + 6.400000e-02, 3.462127e-04, 2.424721e-02, 2.486736e-02, 4.876836e-02, 9.922374e-03, -9.922372e-03 + 6.500000e-02, 3.377481e-04, 2.422993e-02, 2.488517e-02, 4.877734e-02, 9.927796e-03, -9.927795e-03 + 6.600000e-02, 3.287944e-04, 2.421395e-02, 2.490176e-02, 4.878692e-02, 9.923402e-03, -9.923402e-03 + 6.700000e-02, 3.193870e-04, 2.419935e-02, 2.491708e-02, 4.879704e-02, 9.909197e-03, -9.909197e-03 + 6.800000e-02, 3.095644e-04, 2.418618e-02, 2.493106e-02, 4.880768e-02, 9.895821e-03, -9.895810e-03 + 6.900000e-02, 2.993681e-04, 2.417451e-02, 2.494364e-02, 4.881878e-02, 9.911710e-03, -9.911688e-03 + 7.000000e-02, 2.888420e-04, 2.416437e-02, 2.495476e-02, 4.883029e-02, 9.917797e-03, -9.917767e-03 + 7.100000e-02, 2.780313e-04, 2.415582e-02, 2.496436e-02, 4.884215e-02, 9.914069e-03, -9.914040e-03 + 7.200000e-02, 2.669813e-04, 2.414890e-02, 2.497240e-02, 4.885431e-02, 9.900526e-03, -9.900510e-03 + 7.300000e-02, 2.557375e-04, 2.414363e-02, 2.497883e-02, 4.886672e-02, 9.891787e-03, -9.891787e-03 + 7.400000e-02, 2.443447e-04, 2.414005e-02, 2.498362e-02, 4.887933e-02, 9.908385e-03, -9.908385e-03 + 7.500000e-02, 2.328481e-04, 2.413817e-02, 2.498676e-02, 4.889208e-02, 9.915238e-03, -9.915185e-03 + 7.600000e-02, 2.212935e-04, 2.413801e-02, 2.498821e-02, 4.890493e-02, 9.912353e-03, -9.912183e-03 + 7.700000e-02, 2.097280e-04, 2.413958e-02, 2.498799e-02, 4.891784e-02, 9.899664e-03, -9.899379e-03 + 7.800000e-02, 1.981994e-04, 2.414287e-02, 2.498607e-02, 4.893073e-02, 9.895474e-03, -9.895474e-03 + 7.900000e-02, 1.867562e-04, 2.414788e-02, 2.498245e-02, 4.894357e-02, 9.912790e-03, -9.912790e-03 + 8.000000e-02, 1.754469e-04, 2.415459e-02, 2.497714e-02, 4.895628e-02, 9.920304e-03, -9.920304e-03 + 8.100000e-02, 1.643188e-04, 2.416298e-02, 2.497014e-02, 4.896881e-02, 9.918010e-03, -9.918010e-03 + 8.200000e-02, 1.534181e-04, 2.417304e-02, 2.496148e-02, 4.898110e-02, 9.905909e-03, -9.905909e-03 + 8.300000e-02, 1.427892e-04, 2.418472e-02, 2.495118e-02, 4.899311e-02, 9.909187e-03, -9.906520e-03 + 8.400000e-02, 1.324752e-04, 2.419798e-02, 2.493927e-02, 4.900478e-02, 9.926714e-03, -9.924477e-03 + 8.500000e-02, 1.225179e-04, 2.421279e-02, 2.492580e-02, 4.901606e-02, 9.933772e-03, -9.932622e-03 + 8.600000e-02, 1.129582e-04, 2.422908e-02, 2.491081e-02, 4.902692e-02, 9.930946e-03, -9.930946e-03 + 8.700000e-02, 1.038360e-04, 2.424679e-02, 2.489435e-02, 4.903731e-02, 9.919450e-03, -9.919451e-03 + 8.800000e-02, 9.518985e-05, 2.426588e-02, 2.487649e-02, 4.904717e-02, 9.923852e-03, -9.923852e-03 + 8.900000e-02, 8.705675e-05, 2.428625e-02, 2.485727e-02, 4.905647e-02, 9.942511e-03, -9.942308e-03 + 9.000000e-02, 7.947142e-05, 2.430785e-02, 2.483678e-02, 4.906515e-02, 9.954406e-03, -9.950933e-03 + 9.100000e-02, 7.246611e-05, 2.433058e-02, 2.481507e-02, 4.907319e-02, 9.956284e-03, -9.949720e-03 + 9.200000e-02, 6.607044e-05, 2.435437e-02, 2.479224e-02, 4.908054e-02, 9.947641e-03, -9.938669e-03 + 9.300000e-02, 6.031138e-05, 2.437913e-02, 2.476836e-02, 4.908717e-02, 9.948470e-03, -9.945783e-03 + 9.400000e-02, 5.521347e-05, 2.440475e-02, 2.474352e-02, 4.909306e-02, 9.964540e-03, -9.964540e-03 + 9.500000e-02, 5.079878e-05, 2.443115e-02, 2.471781e-02, 4.909817e-02, 9.973445e-03, -9.973445e-03 + 9.600000e-02, 4.708697e-05, 2.445823e-02, 2.469134e-02, 4.910248e-02, 9.972488e-03, -9.972489e-03 + 9.700000e-02, 4.409509e-05, 2.448587e-02, 2.466420e-02, 4.910598e-02, 9.961672e-03, -9.961672e-03 + 9.800000e-02, 4.183737e-05, 2.451398e-02, 2.463649e-02, 4.910863e-02, 9.970172e-03, -9.970172e-03 + 9.900000e-02, 4.032509e-05, 2.454244e-02, 2.460832e-02, 4.911044e-02, 9.989000e-03, -9.989000e-03 + 1.000000e-01, 3.956640e-05, 2.457115e-02, 2.457979e-02, 4.911137e-02, 9.997951e-03, -9.997952e-03 + 1.010000e-01, 3.956640e-05, 2.459999e-02, 2.455102e-02, 4.911143e-02, 9.997017e-03, -1.000000e-02 + 1.020000e-01, 4.032720e-05, 2.462884e-02, 2.452210e-02, 4.911062e-02, 9.986198e-03, -9.995066e-03 + 1.030000e-01, 4.184800e-05, 2.465761e-02, 2.449316e-02, 4.910893e-02, 9.998170e-03, -9.994642e-03 + 1.040000e-01, 4.412514e-05, 2.468617e-02, 2.446431e-02, 4.910635e-02, 1.002060e-02, -1.001330e-02 + 1.050000e-01, 4.715200e-05, 2.471441e-02, 2.443566e-02, 4.910291e-02, 1.003337e-02, -1.002205e-02 + 1.060000e-01, 5.091895e-05, 2.474221e-02, 2.440731e-02, 4.909861e-02, 1.003656e-02, -1.002090e-02 + 1.070000e-01, 5.541328e-05, 2.476948e-02, 2.437939e-02, 4.909345e-02, 1.003468e-02, -1.000984e-02 + 1.080000e-01, 6.061926e-05, 2.479609e-02, 2.435200e-02, 4.908746e-02, 1.006776e-02, -1.001680e-02 + 1.090000e-01, 6.651822e-05, 2.482194e-02, 2.432524e-02, 4.908066e-02, 1.009076e-02, -1.003506e-02 + 1.100000e-01, 7.308883e-05, 2.484692e-02, 2.429923e-02, 4.907307e-02, 1.010377e-02, -1.004339e-02 + 1.110000e-01, 8.030728e-05, 2.487094e-02, 2.427408e-02, 4.906471e-02, 1.010690e-02, -1.004179e-02 + 1.120000e-01, 8.814745e-05, 2.489389e-02, 2.424987e-02, 4.905562e-02, 1.010175e-02, -1.003026e-02 + 1.130000e-01, 9.658092e-05, 2.491568e-02, 2.422672e-02, 4.904582e-02, 1.013393e-02, -1.003453e-02 + 1.140000e-01, 1.055769e-04, 2.493622e-02, 2.420470e-02, 4.903535e-02, 1.015593e-02, -1.005215e-02 + 1.150000e-01, 1.151022e-04, 2.495543e-02, 2.418392e-02, 4.902424e-02, 1.016780e-02, -1.005987e-02 + 1.160000e-01, 1.251212e-04, 2.497321e-02, 2.416444e-02, 4.901254e-02, 1.016964e-02, -1.005764e-02 + 1.170000e-01, 1.355961e-04, 2.498950e-02, 2.414637e-02, 4.900028e-02, 1.016153e-02, -1.004546e-02 + 1.180000e-01, 1.464869e-04, 2.500423e-02, 2.412977e-02, 4.898751e-02, 1.018582e-02, -1.004633e-02 + 1.190000e-01, 1.577523e-04, 2.501733e-02, 2.411471e-02, 4.897429e-02, 1.020606e-02, -1.006337e-02 + 1.200000e-01, 1.693498e-04, 2.502874e-02, 2.410126e-02, 4.896065e-02, 1.021610e-02, -1.007036e-02 + 1.210000e-01, 1.812357e-04, 2.503841e-02, 2.408948e-02, 4.894666e-02, 1.021601e-02, -1.006726e-02 + 1.220000e-01, 1.933658e-04, 2.504630e-02, 2.407942e-02, 4.893235e-02, 1.020588e-02, -1.005406e-02 + 1.230000e-01, 2.056946e-04, 2.505236e-02, 2.407112e-02, 4.891778e-02, 1.021748e-02, -1.005018e-02 + 1.240000e-01, 2.181759e-04, 2.505656e-02, 2.406462e-02, 4.890300e-02, 1.023546e-02, -1.006628e-02 + 1.250000e-01, 2.307622e-04, 2.505889e-02, 2.405995e-02, 4.888807e-02, 1.024321e-02, -1.007267e-02 + 1.260000e-01, 2.434053e-04, 2.505932e-02, 2.405713e-02, 4.887304e-02, 1.024078e-02, -1.006979e-02 + 1.270000e-01, 2.560565e-04, 2.505783e-02, 2.405619e-02, 4.885797e-02, 1.022825e-02, -1.005693e-02 + 1.280000e-01, 2.686669e-04, 2.505444e-02, 2.405714e-02, 4.884292e-02, 1.022525e-02, -1.004666e-02 + 1.290000e-01, 2.811885e-04, 2.504914e-02, 2.405998e-02, 4.882794e-02, 1.024071e-02, -1.006195e-02 + 1.300000e-01, 2.935736e-04, 2.504195e-02, 2.406471e-02, 4.881309e-02, 1.024593e-02, -1.006730e-02 + 1.310000e-01, 3.057757e-04, 2.503289e-02, 2.407131e-02, 4.879842e-02, 1.024096e-02, -1.006269e-02 + 1.320000e-01, 3.177492e-04, 2.502197e-02, 2.407977e-02, 4.878399e-02, 1.022586e-02, -1.004813e-02 + 1.330000e-01, 3.294492e-04, 2.500925e-02, 2.409005e-02, 4.876985e-02, 1.020779e-02, -1.004220e-02 + 1.340000e-01, 3.408311e-04, 2.499475e-02, 2.410212e-02, 4.875605e-02, 1.022072e-02, -1.005812e-02 + 1.350000e-01, 3.518515e-04, 2.497854e-02, 2.411595e-02, 4.874263e-02, 1.022343e-02, -1.006300e-02 + 1.360000e-01, 3.624679e-04, 2.496065e-02, 2.413148e-02, 4.872966e-02, 1.021597e-02, -1.005680e-02 + 1.370000e-01, 3.726391e-04, 2.494117e-02, 2.414865e-02, 4.871718e-02, 1.019840e-02, -1.003971e-02 + 1.380000e-01, 3.823261e-04, 2.492014e-02, 2.416742e-02, 4.870524e-02, 1.017079e-02, -1.001862e-02 + 1.390000e-01, 3.914924e-04, 2.489766e-02, 2.418771e-02, 4.869389e-02, 1.017694e-02, -1.003259e-02 + 1.400000e-01, 4.001038e-04, 2.487380e-02, 2.420946e-02, 4.868316e-02, 1.017742e-02, -1.003664e-02 + 1.410000e-01, 4.081286e-04, 2.484865e-02, 2.423257e-02, 4.867309e-02, 1.016777e-02, -1.003076e-02 + 1.420000e-01, 4.155375e-04, 2.482231e-02, 2.425695e-02, 4.866372e-02, 1.014805e-02, -1.001497e-02 + 1.430000e-01, 4.223031e-04, 2.479486e-02, 2.428253e-02, 4.865508e-02, 1.011832e-02, -9.998188e-03 + 1.440000e-01, 4.283999e-04, 2.476641e-02, 2.430919e-02, 4.864720e-02, 1.011324e-02, -1.001483e-02 + 1.450000e-01, 4.338050e-04, 2.473707e-02, 2.433685e-02, 4.864011e-02, 1.011199e-02, -1.002171e-02 + 1.460000e-01, 4.384975e-04, 2.470695e-02, 2.436538e-02, 4.863383e-02, 1.010068e-02, -1.001857e-02 + 1.470000e-01, 4.424600e-04, 2.467616e-02, 2.439470e-02, 4.862840e-02, 1.007935e-02, -1.000515e-02 + 1.480000e-01, 4.456781e-04, 2.464481e-02, 2.442469e-02, 4.862383e-02, 1.004808e-02, -9.990924e-03 + 1.490000e-01, 4.481408e-04, 2.461304e-02, 2.445523e-02, 4.862013e-02, 1.003557e-02, -1.000376e-02 + 1.500000e-01, 4.498406e-04, 2.458096e-02, 2.448620e-02, 4.861732e-02, 1.003328e-02, -1.000649e-02 + 1.510000e-01, 4.507731e-04, 2.454869e-02, 2.451748e-02, 4.861540e-02, 1.002100e-02, -9.999048e-03 + 1.520000e-01, 4.509364e-04, 2.451637e-02, 2.454894e-02, 4.861437e-02, 9.998796e-03, -9.981324e-03 + 1.530000e-01, 4.503314e-04, 2.448410e-02, 2.458047e-02, 4.861424e-02, 9.966719e-03, -9.953223e-03 + 1.540000e-01, 4.489617e-04, 2.445203e-02, 2.461193e-02, 4.861500e-02, 9.961931e-03, -9.961904e-03 + 1.550000e-01, 4.468333e-04, 2.442027e-02, 2.464321e-02, 4.861664e-02, 9.965606e-03, -9.965546e-03 + 1.560000e-01, 4.439555e-04, 2.438895e-02, 2.467418e-02, 4.861917e-02, 9.959422e-03, -9.959341e-03 + 1.570000e-01, 4.403408e-04, 2.435820e-02, 2.470471e-02, 4.862257e-02, 9.943373e-03, -9.943297e-03 + 1.580000e-01, 4.360052e-04, 2.432813e-02, 2.473470e-02, 4.862682e-02, 9.926465e-03, -9.926465e-03 + 1.590000e-01, 4.309677e-04, 2.429887e-02, 2.476401e-02, 4.863191e-02, 9.940267e-03, -9.940267e-03 + 1.600000e-01, 4.252504e-04, 2.427053e-02, 2.479253e-02, 4.863781e-02, 9.944241e-03, -9.944244e-03 + 1.610000e-01, 4.188779e-04, 2.424322e-02, 2.482014e-02, 4.864448e-02, 9.938382e-03, -9.938400e-03 + 1.620000e-01, 4.118772e-04, 2.421706e-02, 2.484672e-02, 4.865190e-02, 9.922779e-03, -9.922731e-03 + 1.630000e-01, 4.042771e-04, 2.419216e-02, 2.487216e-02, 4.866004e-02, 9.909308e-03, -9.908997e-03 + 1.640000e-01, 3.961088e-04, 2.416861e-02, 2.489636e-02, 4.866886e-02, 9.923304e-03, -9.923304e-03 + 1.650000e-01, 3.874054e-04, 2.414651e-02, 2.491923e-02, 4.867833e-02, 9.927800e-03, -9.927815e-03 + 1.660000e-01, 3.782025e-04, 2.412595e-02, 2.494066e-02, 4.868841e-02, 9.922479e-03, -9.922513e-03 + 1.670000e-01, 3.685380e-04, 2.410702e-02, 2.496057e-02, 4.869905e-02, 9.907348e-03, -9.907397e-03 + 1.680000e-01, 3.584521e-04, 2.408980e-02, 2.497887e-02, 4.871021e-02, 9.899312e-03, -9.897729e-03 + 1.690000e-01, 3.479866e-04, 2.407436e-02, 2.499548e-02, 4.872185e-02, 9.915057e-03, -9.912693e-03 + 1.700000e-01, 3.371851e-04, 2.406076e-02, 2.501032e-02, 4.873390e-02, 9.920626e-03, -9.917866e-03 + 1.710000e-01, 3.260923e-04, 2.404908e-02, 2.502333e-02, 4.874632e-02, 9.915846e-03, -9.913231e-03 + 1.720000e-01, 3.147534e-04, 2.403935e-02, 2.503446e-02, 4.875906e-02, 9.900649e-03, -9.898786e-03 + 1.730000e-01, 3.032148e-04, 2.403163e-02, 2.504364e-02, 4.877205e-02, 9.893774e-03, -9.893780e-03 + 1.740000e-01, 2.915232e-04, 2.402596e-02, 2.505083e-02, 4.878526e-02, 9.909447e-03, -9.909475e-03 + 1.750000e-01, 2.797262e-04, 2.402235e-02, 2.505600e-02, 4.879862e-02, 9.915321e-03, -9.915373e-03 + 1.760000e-01, 2.678720e-04, 2.402084e-02, 2.505912e-02, 4.881208e-02, 9.911392e-03, -9.911460e-03 + 1.770000e-01, 2.560094e-04, 2.402143e-02, 2.506016e-02, 4.882559e-02, 9.899148e-03, -9.897731e-03 + 1.780000e-01, 2.441875e-04, 2.402414e-02, 2.505912e-02, 4.883907e-02, 9.902387e-03, -9.897553e-03 + 1.790000e-01, 2.324552e-04, 2.402895e-02, 2.505599e-02, 4.885248e-02, 9.916595e-03, -9.913970e-03 + 1.800000e-01, 2.208608e-04, 2.403586e-02, 2.505076e-02, 4.886577e-02, 9.920665e-03, -9.920579e-03 + 1.810000e-01, 2.094523e-04, 2.404485e-02, 2.504346e-02, 4.887886e-02, 9.917286e-03, -9.917365e-03 + 1.820000e-01, 1.982763e-04, 2.405589e-02, 2.503409e-02, 4.889171e-02, 9.904259e-03, -9.904323e-03 + 1.830000e-01, 1.873787e-04, 2.406894e-02, 2.502270e-02, 4.890426e-02, 9.908636e-03, -9.908684e-03 + 1.840000e-01, 1.768041e-04, 2.408396e-02, 2.500930e-02, 4.891646e-02, 9.925667e-03, -9.925742e-03 + 1.850000e-01, 1.665962e-04, 2.410090e-02, 2.499395e-02, 4.892826e-02, 9.932884e-03, -9.932973e-03 + 1.860000e-01, 1.567973e-04, 2.411969e-02, 2.497671e-02, 4.893960e-02, 9.930280e-03, -9.930361e-03 + 1.870000e-01, 1.474483e-04, 2.414028e-02, 2.495762e-02, 4.895045e-02, 9.918424e-03, -9.917902e-03 + 1.880000e-01, 1.385885e-04, 2.416258e-02, 2.493676e-02, 4.896075e-02, 9.936195e-03, -9.926092e-03 + 1.890000e-01, 1.302551e-04, 2.418651e-02, 2.491419e-02, 4.897045e-02, 9.954142e-03, -9.943642e-03 + 1.900000e-01, 1.224831e-04, 2.421199e-02, 2.489001e-02, 4.897952e-02, 9.961486e-03, -9.951338e-03 + 1.910000e-01, 1.153053e-04, 2.423892e-02, 2.486428e-02, 4.898790e-02, 9.958231e-03, -9.949165e-03 + 1.920000e-01, 1.087519e-04, 2.426721e-02, 2.483712e-02, 4.899558e-02, 9.944498e-03, -9.937126e-03 + 1.930000e-01, 1.028508e-04, 2.429673e-02, 2.480862e-02, 4.900250e-02, 9.947976e-03, -9.948080e-03 + 1.940000e-01, 9.762740e-05, 2.432739e-02, 2.477889e-02, 4.900865e-02, 9.965803e-03, -9.965917e-03 + 1.950000e-01, 9.310448e-05, 2.435906e-02, 2.474803e-02, 4.901399e-02, 9.973776e-03, -9.973869e-03 + 1.960000e-01, 8.930214e-05, 2.439163e-02, 2.471617e-02, 4.901849e-02, 9.971888e-03, -9.971924e-03 + 1.970000e-01, 8.623765e-05, 2.442496e-02, 2.468342e-02, 4.902214e-02, 9.960140e-03, -9.960140e-03 + 1.980000e-01, 8.392535e-05, 2.445894e-02, 2.464990e-02, 4.902492e-02, 9.972374e-03, -9.972500e-03 + 1.990000e-01, 8.237658e-05, 2.449342e-02, 2.461576e-02, 4.902680e-02, 9.990269e-03, -9.990385e-03 + 2.000000e-01, 8.159958e-05, 2.452828e-02, 2.458110e-02, 4.902778e-02, 9.998286e-03, -9.998353e-03 + 2.010000e-01, 8.159958e-05, 2.456338e-02, 2.454607e-02, 4.902786e-02, 1.000000e-02, -9.996417e-03 + 2.020000e-01, 8.237876e-05, 2.459858e-02, 2.451081e-02, 4.902701e-02, 9.995066e-03, -9.984665e-03 + 2.030000e-01, 8.393626e-05, 2.463374e-02, 2.447545e-02, 4.902526e-02, 1.000235e-02, -1.000176e-02 + 2.040000e-01, 8.626820e-05, 2.466873e-02, 2.444012e-02, 4.902258e-02, 1.002425e-02, -1.002561e-02 + 2.050000e-01, 8.936758e-05, 2.470340e-02, 2.440498e-02, 4.901900e-02, 1.003656e-02, -1.003950e-02 + 2.060000e-01, 9.322436e-05, 2.473761e-02, 2.437014e-02, 4.901453e-02, 1.003915e-02, -1.004344e-02 + 2.070000e-01, 9.782541e-05, 2.477124e-02, 2.433575e-02, 4.900917e-02, 1.003833e-02, -1.005148e-02 + 2.080000e-01, 1.031547e-04, 2.480414e-02, 2.430195e-02, 4.900294e-02, 1.007028e-02, -1.008474e-02 + 2.090000e-01, 1.091931e-04, 2.483618e-02, 2.426887e-02, 4.899586e-02, 1.009250e-02, -1.010799e-02 + 2.100000e-01, 1.159192e-04, 2.486724e-02, 2.423665e-02, 4.898797e-02, 1.010496e-02, -1.012119e-02 + 2.110000e-01, 1.233084e-04, 2.489718e-02, 2.420540e-02, 4.897928e-02, 1.010758e-02, -1.012436e-02 + 2.120000e-01, 1.313340e-04, 2.492588e-02, 2.417527e-02, 4.896982e-02, 1.010492e-02, -1.012646e-02 + 2.130000e-01, 1.399665e-04, 2.495323e-02, 2.414636e-02, 4.895962e-02, 1.013588e-02, -1.015838e-02 + 2.140000e-01, 1.491739e-04, 2.497911e-02, 2.411879e-02, 4.894873e-02, 1.015696e-02, -1.018021e-02 + 2.150000e-01, 1.589220e-04, 2.500342e-02, 2.409269e-02, 4.893718e-02, 1.016817e-02, -1.019192e-02 + 2.160000e-01, 1.691742e-04, 2.502604e-02, 2.406815e-02, 4.892502e-02, 1.016948e-02, -1.019353e-02 + 2.170000e-01, 1.798918e-04, 2.504689e-02, 2.404528e-02, 4.891228e-02, 1.016081e-02, -1.018523e-02 + 2.180000e-01, 1.910345e-04, 2.506588e-02, 2.402417e-02, 4.889901e-02, 1.018772e-02, -1.021509e-02 + 2.190000e-01, 2.025602e-04, 2.508292e-02, 2.400491e-02, 4.888527e-02, 1.020691e-02, -1.023482e-02 + 2.200000e-01, 2.144255e-04, 2.509793e-02, 2.398759e-02, 4.887110e-02, 1.021613e-02, -1.024438e-02 + 2.210000e-01, 2.265858e-04, 2.511086e-02, 2.397228e-02, 4.885656e-02, 1.021541e-02, -1.024378e-02 + 2.220000e-01, 2.389952e-04, 2.512164e-02, 2.395904e-02, 4.884169e-02, 1.020469e-02, -1.023306e-02 + 2.230000e-01, 2.516069e-04, 2.513022e-02, 2.394794e-02, 4.882655e-02, 1.021935e-02, -1.024886e-02 + 2.240000e-01, 2.643728e-04, 2.513655e-02, 2.393902e-02, 4.881120e-02, 1.023621e-02, -1.026598e-02 + 2.250000e-01, 2.772444e-04, 2.514061e-02, 2.393233e-02, 4.879570e-02, 1.024302e-02, -1.027292e-02 + 2.260000e-01, 2.901723e-04, 2.514236e-02, 2.392791e-02, 4.878010e-02, 1.023984e-02, -1.026966e-02 + 2.270000e-01, 3.031071e-04, 2.514180e-02, 2.392577e-02, 4.876446e-02, 1.022666e-02, -1.025626e-02 + 2.280000e-01, 3.159995e-04, 2.513890e-02, 2.392593e-02, 4.874884e-02, 1.022703e-02, -1.025582e-02 + 2.290000e-01, 3.288004e-04, 2.513369e-02, 2.392841e-02, 4.873330e-02, 1.024133e-02, -1.027010e-02 + 2.300000e-01, 3.414614e-04, 2.512615e-02, 2.393320e-02, 4.871789e-02, 1.024555e-02, -1.027418e-02 + 2.310000e-01, 3.539344e-04, 2.511633e-02, 2.394029e-02, 4.870268e-02, 1.023974e-02, -1.026807e-02 + 2.320000e-01, 3.661722e-04, 2.510423e-02, 2.394965e-02, 4.868771e-02, 1.022392e-02, -1.025181e-02 + 2.330000e-01, 3.781282e-04, 2.508992e-02, 2.396126e-02, 4.867305e-02, 1.020944e-02, -1.023472e-02 + 2.340000e-01, 3.897568e-04, 2.507342e-02, 2.397509e-02, 4.865874e-02, 1.022122e-02, -1.024619e-02 + 2.350000e-01, 4.010133e-04, 2.505479e-02, 2.399107e-02, 4.864485e-02, 1.022289e-02, -1.024747e-02 + 2.360000e-01, 4.118547e-04, 2.503411e-02, 2.400917e-02, 4.863142e-02, 1.021452e-02, -1.023860e-02 + 2.370000e-01, 4.222395e-04, 2.501144e-02, 2.402931e-02, 4.861851e-02, 1.019615e-02, -1.021959e-02 + 2.380000e-01, 4.321284e-04, 2.498686e-02, 2.405143e-02, 4.860616e-02, 1.016779e-02, -1.019050e-02 + 2.390000e-01, 4.414842e-04, 2.496046e-02, 2.407544e-02, 4.859442e-02, 1.017731e-02, -1.019608e-02 + 2.400000e-01, 4.502717e-04, 2.493234e-02, 2.410126e-02, 4.858333e-02, 1.017673e-02, -1.019490e-02 + 2.410000e-01, 4.584583e-04, 2.490260e-02, 2.412879e-02, 4.857293e-02, 1.016613e-02, -1.018362e-02 + 2.420000e-01, 4.660131e-04, 2.487135e-02, 2.415792e-02, 4.856326e-02, 1.014556e-02, -1.016226e-02 + 2.430000e-01, 4.729079e-04, 2.483871e-02, 2.418854e-02, 4.855435e-02, 1.011504e-02, -1.013086e-02 + 2.440000e-01, 4.791165e-04, 2.480480e-02, 2.422055e-02, 4.854623e-02, 1.011350e-02, -1.012425e-02 + 2.450000e-01, 4.846157e-04, 2.476975e-02, 2.425381e-02, 4.853894e-02, 1.011120e-02, -1.012120e-02 + 2.460000e-01, 4.893850e-04, 2.473368e-02, 2.428821e-02, 4.853251e-02, 1.009891e-02, -1.010811e-02 + 2.470000e-01, 4.934068e-04, 2.469675e-02, 2.432361e-02, 4.852695e-02, 1.007670e-02, -1.008502e-02 + 2.480000e-01, 4.966671e-04, 2.465908e-02, 2.435987e-02, 4.852228e-02, 1.004460e-02, -1.005196e-02 + 2.490000e-01, 4.991546e-04, 2.462082e-02, 2.439686e-02, 4.851853e-02, 1.003576e-02, -1.003747e-02 + 2.500000e-01, 5.008613e-04, 2.458212e-02, 2.443443e-02, 4.851569e-02, 1.003243e-02, -1.003331e-02 + 2.510000e-01, 5.017823e-04, 2.454314e-02, 2.447243e-02, 4.851378e-02, 1.001917e-02, -1.001921e-02 + 2.520000e-01, 5.019152e-04, 2.450401e-02, 2.451071e-02, 4.851280e-02, 9.996053e-03, -9.995198e-03 + 2.530000e-01, 5.012609e-04, 2.446490e-02, 2.454913e-02, 4.851276e-02, 9.963118e-03, -9.961305e-03 + 2.540000e-01, 4.998229e-04, 2.442595e-02, 2.458752e-02, 4.851365e-02, 9.962902e-03, -9.962669e-03 + 2.550000e-01, 4.976082e-04, 2.438732e-02, 2.462575e-02, 4.851546e-02, 9.966608e-03, -9.965381e-03 + 2.560000e-01, 4.946267e-04, 2.434916e-02, 2.466365e-02, 4.851819e-02, 9.960365e-03, -9.958443e-03 + 2.570000e-01, 4.908919e-04, 2.431163e-02, 2.470109e-02, 4.852182e-02, 9.943959e-03, -9.942059e-03 + 2.580000e-01, 4.864201e-04, 2.427487e-02, 2.473789e-02, 4.852634e-02, 9.928190e-03, -9.928190e-03 + 2.590000e-01, 4.812309e-04, 2.423902e-02, 2.477393e-02, 4.853173e-02, 9.941063e-03, -9.941063e-03 + 2.600000e-01, 4.753464e-04, 2.420424e-02, 2.480905e-02, 4.853795e-02, 9.944108e-03, -9.944108e-03 + 2.610000e-01, 4.687915e-04, 2.417066e-02, 2.484311e-02, 4.854498e-02, 9.937320e-03, -9.937927e-03 + 2.620000e-01, 4.615935e-04, 2.413842e-02, 2.487596e-02, 4.855279e-02, 9.920707e-03, -9.921933e-03 + 2.630000e-01, 4.537819e-04, 2.410765e-02, 2.490748e-02, 4.856134e-02, 9.916312e-03, -9.910769e-03 + 2.640000e-01, 4.453889e-04, 2.407847e-02, 2.493753e-02, 4.857061e-02, 9.929424e-03, -9.924150e-03 + 2.650000e-01, 4.364490e-04, 2.405100e-02, 2.496600e-02, 4.858055e-02, 9.931842e-03, -9.928090e-03 + 2.660000e-01, 4.269991e-04, 2.402536e-02, 2.499275e-02, 4.859112e-02, 9.923775e-03, -9.922522e-03 + 2.670000e-01, 4.170782e-04, 2.400166e-02, 2.501769e-02, 4.860227e-02, 9.905558e-03, -9.907036e-03 + 2.680000e-01, 4.067274e-04, 2.397999e-02, 2.504070e-02, 4.861396e-02, 9.899561e-03, -9.899560e-03 + 2.690000e-01, 3.959893e-04, 2.396045e-02, 2.506168e-02, 4.862614e-02, 9.913589e-03, -9.913702e-03 + 2.700000e-01, 3.849083e-04, 2.394312e-02, 2.508054e-02, 4.863875e-02, 9.917814e-03, -9.918648e-03 + 2.710000e-01, 3.735296e-04, 2.392807e-02, 2.509720e-02, 4.865174e-02, 9.912233e-03, -9.913737e-03 + 2.720000e-01, 3.618996e-04, 2.391538e-02, 2.511158e-02, 4.866506e-02, 9.897084e-03, -9.898808e-03 + 2.730000e-01, 3.500658e-04, 2.390509e-02, 2.512362e-02, 4.867865e-02, 9.900708e-03, -9.895663e-03 + 2.740000e-01, 3.380764e-04, 2.389727e-02, 2.513326e-02, 4.869245e-02, 9.915680e-03, -9.910996e-03 + 2.750000e-01, 3.259805e-04, 2.389194e-02, 2.514045e-02, 4.870641e-02, 9.920499e-03, -9.916695e-03 + 2.760000e-01, 3.138276e-04, 2.388915e-02, 2.514515e-02, 4.872047e-02, 9.915053e-03, -9.912427e-03 + 2.770000e-01, 3.016678e-04, 2.388890e-02, 2.514735e-02, 4.873457e-02, 9.899289e-03, -9.898034e-03 + 2.780000e-01, 2.895510e-04, 2.389120e-02, 2.514700e-02, 4.874865e-02, 9.899497e-03, -9.899782e-03 + 2.790000e-01, 2.775271e-04, 2.389607e-02, 2.514411e-02, 4.876265e-02, 9.914961e-03, -9.916054e-03 + 2.800000e-01, 2.656454e-04, 2.390347e-02, 2.513868e-02, 4.877651e-02, 9.920622e-03, -9.922424e-03 + 2.810000e-01, 2.539547e-04, 2.391341e-02, 2.513072e-02, 4.879017e-02, 9.916474e-03, -9.918705e-03 + 2.820000e-01, 2.425028e-04, 2.392584e-02, 2.512024e-02, 4.880358e-02, 9.902521e-03, -9.904756e-03 + 2.830000e-01, 2.313367e-04, 2.394072e-02, 2.510729e-02, 4.881667e-02, 9.912617e-03, -9.911477e-03 + 2.840000e-01, 2.205024e-04, 2.395801e-02, 2.509189e-02, 4.882940e-02, 9.930317e-03, -9.928392e-03 + 2.850000e-01, 2.100446e-04, 2.397764e-02, 2.507411e-02, 4.884171e-02, 9.938262e-03, -9.935270e-03 + 2.860000e-01, 2.000068e-04, 2.399955e-02, 2.505400e-02, 4.885354e-02, 9.936398e-03, -9.931933e-03 + 2.870000e-01, 1.904308e-04, 2.402365e-02, 2.503163e-02, 4.886485e-02, 9.924618e-03, -9.918280e-03 + 2.880000e-01, 1.813564e-04, 2.404986e-02, 2.500708e-02, 4.887558e-02, 9.936993e-03, -9.929473e-03 + 2.890000e-01, 1.728216e-04, 2.407808e-02, 2.498043e-02, 4.888569e-02, 9.951932e-03, -9.946807e-03 + 2.900000e-01, 1.648620e-04, 2.410821e-02, 2.495179e-02, 4.889514e-02, 9.956771e-03, -9.953960e-03 + 2.910000e-01, 1.575112e-04, 2.414014e-02, 2.492126e-02, 4.890389e-02, 9.951654e-03, -9.950787e-03 + 2.920000e-01, 1.508001e-04, 2.417374e-02, 2.488894e-02, 4.891189e-02, 9.936685e-03, -9.937240e-03 + 2.930000e-01, 1.447574e-04, 2.420889e-02, 2.485497e-02, 4.891910e-02, 9.950082e-03, -9.952022e-03 + 2.940000e-01, 1.394089e-04, 2.424546e-02, 2.481946e-02, 4.892551e-02, 9.966978e-03, -9.969491e-03 + 2.950000e-01, 1.347780e-04, 2.428330e-02, 2.478255e-02, 4.893107e-02, 9.974020e-03, -9.976645e-03 + 2.960000e-01, 1.308851e-04, 2.432227e-02, 2.474439e-02, 4.893577e-02, 9.971200e-03, -9.973386e-03 + 2.970000e-01, 1.277479e-04, 2.436222e-02, 2.470510e-02, 4.893957e-02, 9.958520e-03, -9.959738e-03 + 2.980000e-01, 1.253809e-04, 2.440300e-02, 2.466485e-02, 4.894247e-02, 9.974487e-03, -9.976922e-03 + 2.990000e-01, 1.237955e-04, 2.444444e-02, 2.462379e-02, 4.894443e-02, 9.991449e-03, -9.994212e-03 + 3.000000e-01, 1.230001e-04, 2.448639e-02, 2.458207e-02, 4.894546e-02, 9.998532e-03, -1.000107e-02 + 3.010000e-01, 1.230001e-04, 2.452869e-02, 2.453986e-02, 4.894554e-02, 9.995730e-03, -1.000000e-02 + 3.020000e-01, 1.237977e-04, 2.457116e-02, 2.449732e-02, 4.894467e-02, 9.983044e-03, -9.995066e-03 + 3.030000e-01, 1.253918e-04, 2.461363e-02, 2.445461e-02, 4.894285e-02, 1.001093e-02, -1.000441e-02 + 3.040000e-01, 1.277783e-04, 2.465595e-02, 2.441191e-02, 4.894008e-02, 1.003475e-02, -1.002720e-02 + 3.050000e-01, 1.309501e-04, 2.469794e-02, 2.436937e-02, 4.893636e-02, 1.004862e-02, -1.004010e-02 + 3.060000e-01, 1.348966e-04, 2.473944e-02, 2.432718e-02, 4.893171e-02, 1.005256e-02, -1.004314e-02 + 3.070000e-01, 1.396046e-04, 2.478027e-02, 2.428548e-02, 4.892615e-02, 1.006937e-02, -1.005519e-02 + 3.080000e-01, 1.450573e-04, 2.482028e-02, 2.424446e-02, 4.891968e-02, 1.010246e-02, -1.008740e-02 + 3.090000e-01, 1.512356e-04, 2.485930e-02, 2.420427e-02, 4.891233e-02, 1.012554e-02, -1.010958e-02 + 3.100000e-01, 1.581171e-04, 2.489718e-02, 2.416507e-02, 4.890413e-02, 1.013859e-02, -1.012177e-02 + 3.110000e-01, 1.656768e-04, 2.493376e-02, 2.412702e-02, 4.889510e-02, 1.014162e-02, -1.012401e-02 + 3.120000e-01, 1.738872e-04, 2.496889e-02, 2.409028e-02, 4.888528e-02, 1.015152e-02, -1.013006e-02 + 3.130000e-01, 1.827178e-04, 2.500243e-02, 2.405499e-02, 4.887470e-02, 1.018307e-02, -1.016094e-02 + 3.140000e-01, 1.921360e-04, 2.503424e-02, 2.402129e-02, 4.886339e-02, 1.020453e-02, -1.018172e-02 + 3.150000e-01, 2.021064e-04, 2.506418e-02, 2.398932e-02, 4.885140e-02, 1.021589e-02, -1.019242e-02 + 3.160000e-01, 2.125916e-04, 2.509214e-02, 2.395922e-02, 4.883877e-02, 1.021715e-02, -1.019307e-02 + 3.170000e-01, 2.235523e-04, 2.511799e-02, 2.393111e-02, 4.882555e-02, 1.021526e-02, -1.018864e-02 + 3.180000e-01, 2.349470e-04, 2.514162e-02, 2.390511e-02, 4.881178e-02, 1.024452e-02, -1.021748e-02 + 3.190000e-01, 2.467330e-04, 2.516293e-02, 2.388132e-02, 4.879752e-02, 1.026363e-02, -1.023617e-02 + 3.200000e-01, 2.588656e-04, 2.518183e-02, 2.385986e-02, 4.878282e-02, 1.027256e-02, -1.024471e-02 + 3.210000e-01, 2.712992e-04, 2.519823e-02, 2.384079e-02, 4.876773e-02, 1.027134e-02, -1.024314e-02 + 3.220000e-01, 2.839865e-04, 2.521206e-02, 2.382423e-02, 4.875230e-02, 1.025998e-02, -1.023149e-02 + 3.230000e-01, 2.968795e-04, 2.522326e-02, 2.381022e-02, 4.873660e-02, 1.028031e-02, -1.025102e-02 + 3.240000e-01, 3.099289e-04, 2.523177e-02, 2.379885e-02, 4.872068e-02, 1.029653e-02, -1.026712e-02 + 3.250000e-01, 3.230850e-04, 2.523754e-02, 2.379015e-02, 4.870461e-02, 1.030255e-02, -1.027303e-02 + 3.260000e-01, 3.362974e-04, 2.524055e-02, 2.378419e-02, 4.868844e-02, 1.029838e-02, -1.026879e-02 + 3.270000e-01, 3.495158e-04, 2.524077e-02, 2.378098e-02, 4.867223e-02, 1.028404e-02, -1.025444e-02 + 3.280000e-01, 3.626899e-04, 2.523818e-02, 2.378055e-02, 4.865604e-02, 1.028639e-02, -1.025771e-02 + 3.290000e-01, 3.757695e-04, 2.523280e-02, 2.378291e-02, 4.863994e-02, 1.029948e-02, -1.027098e-02 + 3.300000e-01, 3.887049e-04, 2.522462e-02, 2.378807e-02, 4.862398e-02, 1.030236e-02, -1.027405e-02 + 3.310000e-01, 4.014469e-04, 2.521367e-02, 2.379600e-02, 4.860822e-02, 1.029504e-02, -1.026695e-02 + 3.320000e-01, 4.139469e-04, 2.519999e-02, 2.380669e-02, 4.859273e-02, 1.027755e-02, -1.024973e-02 + 3.330000e-01, 4.261571e-04, 2.518361e-02, 2.382010e-02, 4.857755e-02, 1.026160e-02, -1.023635e-02 + 3.340000e-01, 4.380308e-04, 2.516459e-02, 2.383619e-02, 4.856275e-02, 1.027160e-02, -1.024682e-02 + 3.350000e-01, 4.495225e-04, 2.514299e-02, 2.385491e-02, 4.854838e-02, 1.027143e-02, -1.024710e-02 + 3.360000e-01, 4.605883e-04, 2.511889e-02, 2.387619e-02, 4.853449e-02, 1.026108e-02, -1.023723e-02 + 3.370000e-01, 4.711863e-04, 2.509238e-02, 2.389995e-02, 4.852115e-02, 1.024059e-02, -1.021726e-02 + 3.380000e-01, 4.812761e-04, 2.506355e-02, 2.392611e-02, 4.850839e-02, 1.020999e-02, -1.018850e-02 + 3.390000e-01, 4.908199e-04, 2.503251e-02, 2.395458e-02, 4.849626e-02, 1.021513e-02, -1.019647e-02 + 3.400000e-01, 4.997815e-04, 2.499936e-02, 2.398524e-02, 4.848481e-02, 1.021228e-02, -1.019430e-02 + 3.410000e-01, 5.081273e-04, 2.496423e-02, 2.401799e-02, 4.847409e-02, 1.019932e-02, -1.018203e-02 + 3.420000e-01, 5.158256e-04, 2.492725e-02, 2.405269e-02, 4.846412e-02, 1.017628e-02, -1.015969e-02 + 3.430000e-01, 5.228475e-04, 2.488856e-02, 2.408923e-02, 4.845494e-02, 1.014318e-02, -1.012735e-02 + 3.440000e-01, 5.291664e-04, 2.484831e-02, 2.412745e-02, 4.844660e-02, 1.013515e-02, -1.012446e-02 + 3.450000e-01, 5.347587e-04, 2.480665e-02, 2.416723e-02, 4.843912e-02, 1.013029e-02, -1.012043e-02 + 3.460000e-01, 5.396038e-04, 2.476373e-02, 2.420840e-02, 4.843252e-02, 1.011540e-02, -1.010635e-02 + 3.470000e-01, 5.436841e-04, 2.471972e-02, 2.425080e-02, 4.842683e-02, 1.009050e-02, -1.008229e-02 + 3.480000e-01, 5.469851e-04, 2.467479e-02, 2.429428e-02, 4.842208e-02, 1.005563e-02, -1.004828e-02 + 3.490000e-01, 5.494954e-04, 2.462911e-02, 2.433865e-02, 4.841827e-02, 1.003923e-02, -1.003756e-02 + 3.500000e-01, 5.512067e-04, 2.458286e-02, 2.438376e-02, 4.841541e-02, 1.003322e-02, -1.003243e-02 + 3.510000e-01, 5.521136e-04, 2.453621e-02, 2.442942e-02, 4.841352e-02, 1.001740e-02, -1.001736e-02 + 3.520000e-01, 5.522138e-04, 2.448936e-02, 2.447545e-02, 4.841260e-02, 9.993537e-03, -9.992371e-03 + 3.530000e-01, 5.515081e-04, 2.444248e-02, 2.452168e-02, 4.841265e-02, 9.959787e-03, -9.957530e-03 + 3.540000e-01, 5.500006e-04, 2.439576e-02, 2.456791e-02, 4.841367e-02, 9.958671e-03, -9.958775e-03 + 3.550000e-01, 5.476985e-04, 2.434937e-02, 2.461397e-02, 4.841565e-02, 9.961034e-03, -9.960236e-03 + 3.560000e-01, 5.446125e-04, 2.430351e-02, 2.465968e-02, 4.841858e-02, 9.954029e-03, -9.952689e-03 + 3.570000e-01, 5.407563e-04, 2.425836e-02, 2.470485e-02, 4.842245e-02, 9.937440e-03, -9.936036e-03 + 3.580000e-01, 5.361468e-04, 2.421409e-02, 2.474929e-02, 4.842724e-02, 9.933243e-03, -9.925643e-03 + 3.590000e-01, 5.308038e-04, 2.417088e-02, 2.479284e-02, 4.843292e-02, 9.945307e-03, -9.936902e-03 + 3.600000e-01, 5.247499e-04, 2.412891e-02, 2.483530e-02, 4.843947e-02, 9.947380e-03, -9.939165e-03 + 3.610000e-01, 5.180103e-04, 2.408835e-02, 2.487653e-02, 4.844686e-02, 9.939455e-03, -9.932390e-03 + 3.620000e-01, 5.106130e-04, 2.404935e-02, 2.491633e-02, 4.845507e-02, 9.921495e-03, -9.916401e-03 + 3.630000e-01, 5.025886e-04, 2.401208e-02, 2.495455e-02, 4.846404e-02, 9.908599e-03, -9.907855e-03 + 3.640000e-01, 4.939700e-04, 2.397669e-02, 2.499104e-02, 4.847376e-02, 9.919929e-03, -9.920098e-03 + 3.650000e-01, 4.847929e-04, 2.394333e-02, 2.502563e-02, 4.848417e-02, 9.921911e-03, -9.923331e-03 + 3.660000e-01, 4.750951e-04, 2.391213e-02, 2.505820e-02, 4.849523e-02, 9.914677e-03, -9.917419e-03 + 3.670000e-01, 4.649166e-04, 2.388322e-02, 2.508859e-02, 4.850690e-02, 9.898274e-03, -9.902117e-03 + 3.680000e-01, 4.542995e-04, 2.385672e-02, 2.511670e-02, 4.851912e-02, 9.904881e-03, -9.896398e-03 + 3.690000e-01, 4.432872e-04, 2.383274e-02, 2.514238e-02, 4.853184e-02, 9.918978e-03, -9.909834e-03 + 3.700000e-01, 4.319248e-04, 2.381140e-02, 2.516553e-02, 4.854501e-02, 9.922878e-03, -9.914228e-03 + 3.710000e-01, 4.202587e-04, 2.379278e-02, 2.518605e-02, 4.855857e-02, 9.916633e-03, -9.909345e-03 + 3.720000e-01, 4.083365e-04, 2.377698e-02, 2.520383e-02, 4.857248e-02, 9.900331e-03, -9.894875e-03 + 3.730000e-01, 3.962068e-04, 2.376404e-02, 2.521882e-02, 4.858666e-02, 9.898881e-03, -9.892301e-03 + 3.740000e-01, 3.839191e-04, 2.375404e-02, 2.523094e-02, 4.860106e-02, 9.911198e-03, -9.906663e-03 + 3.750000e-01, 3.715238e-04, 2.374700e-02, 2.524014e-02, 4.861562e-02, 9.913512e-03, -9.911953e-03 + 3.760000e-01, 3.590717e-04, 2.374299e-02, 2.524637e-02, 4.863029e-02, 9.905947e-03, -9.907937e-03 + 3.770000e-01, 3.466137e-04, 2.374203e-02, 2.524957e-02, 4.864499e-02, 9.888686e-03, -9.894323e-03 + 3.780000e-01, 3.342011e-04, 2.374414e-02, 2.524972e-02, 4.865966e-02, 9.894634e-03, -9.897803e-03 + 3.790000e-01, 3.218846e-04, 2.374933e-02, 2.524681e-02, 4.867425e-02, 9.910362e-03, -9.913023e-03 + 3.800000e-01, 3.097146e-04, 2.375756e-02, 2.524084e-02, 4.868869e-02, 9.916448e-03, -9.918715e-03 + 3.810000e-01, 2.977409e-04, 2.376881e-02, 2.523185e-02, 4.870292e-02, 9.912785e-03, -9.914682e-03 + 3.820000e-01, 2.860126e-04, 2.378303e-02, 2.521987e-02, 4.871688e-02, 9.899334e-03, -9.900775e-03 + 3.830000e-01, 2.745778e-04, 2.380020e-02, 2.520490e-02, 4.873052e-02, 9.920402e-03, -9.909026e-03 + 3.840000e-01, 2.634836e-04, 2.382027e-02, 2.518699e-02, 4.874378e-02, 9.937013e-03, -9.926470e-03 + 3.850000e-01, 2.527759e-04, 2.384317e-02, 2.516619e-02, 4.875659e-02, 9.943511e-03, -9.934242e-03 + 3.860000e-01, 2.424990e-04, 2.386883e-02, 2.514258e-02, 4.876891e-02, 9.939761e-03, -9.931851e-03 + 3.870000e-01, 2.326954e-04, 2.389713e-02, 2.511625e-02, 4.878068e-02, 9.925686e-03, -9.918912e-03 + 3.880000e-01, 2.234060e-04, 2.392795e-02, 2.508730e-02, 4.879185e-02, 9.934796e-03, -9.923327e-03 + 3.890000e-01, 2.146693e-04, 2.396119e-02, 2.505585e-02, 4.880238e-02, 9.948750e-03, -9.939998e-03 + 3.900000e-01, 2.065220e-04, 2.399673e-02, 2.502200e-02, 4.881221e-02, 9.952805e-03, -9.947704e-03 + 3.910000e-01, 1.989980e-04, 2.403444e-02, 2.498586e-02, 4.882131e-02, 9.946908e-03, -9.946173e-03 + 3.920000e-01, 1.921292e-04, 2.407420e-02, 2.494756e-02, 4.882963e-02, 9.931014e-03, -9.935057e-03 + 3.930000e-01, 1.859447e-04, 2.411584e-02, 2.490725e-02, 4.883714e-02, 9.936278e-03, -9.953586e-03 + 3.940000e-01, 1.804712e-04, 2.415920e-02, 2.486508e-02, 4.884381e-02, 9.951823e-03, -9.969742e-03 + 3.950000e-01, 1.757323e-04, 2.420412e-02, 2.482121e-02, 4.884960e-02, 9.959632e-03, -9.975561e-03 + 3.960000e-01, 1.717489e-04, 2.425042e-02, 2.477581e-02, 4.885448e-02, 9.958022e-03, -9.971159e-03 + 3.970000e-01, 1.685389e-04, 2.429793e-02, 2.472906e-02, 4.885844e-02, 9.946975e-03, -9.956755e-03 + 3.980000e-01, 1.661171e-04, 2.434644e-02, 2.468114e-02, 4.886146e-02, 9.958852e-03, -9.972743e-03 + 3.990000e-01, 1.644950e-04, 2.439577e-02, 2.463223e-02, 4.886351e-02, 9.980267e-03, -9.991557e-03 + 4.000000e-01, 1.636813e-04, 2.444573e-02, 2.458253e-02, 4.886458e-02, 9.995066e-03, -1.000028e-02 diff --git a/test/test_model/test_common/test_model_solver/test_model_solver.verified b/test/test_model/test_common/test_model_solver/test_model_solver_mumps.verified similarity index 100% copy from test/test_model/test_common/test_model_solver/test_model_solver.verified copy to test/test_model/test_common/test_model_solver/test_model_solver_mumps.verified diff --git a/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh b/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh index 0f88c3f59..8a2fa42d5 100644 --- a/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh +++ b/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh @@ -1,451 +1,450 @@ /** * @file test_model_solver_my_model.hh * * @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 "aka_iterators.hh" #include "boundary_condition.hh" #include "communicator.hh" #include "data_accessor.hh" #include "dof_manager_default.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "model_solver.hh" #include "periodic_node_synchronizer.hh" #include "solver_vector_default.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { #ifndef __AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH__ #define __AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH__ /** * =\o-----o-----o-> F * | | * |---- L ----| */ class MyModel : public ModelSolver, public BoundaryCondition, public DataAccessor { public: - MyModel(Real F, Mesh & mesh, bool lumped) + MyModel(Real F, Mesh & mesh, bool lumped, const ID & dof_manager_type = "default") : ModelSolver(mesh, ModelType::_model, "model_solver", 0), 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") { this->initBC(*this, displacement, forces); - this->initDOFManager("petsc"); - //this->initDOFManager("default"); + this->initDOFManager(dof_manager_type); this->getDOFManager().registerDOFs("disp", displacement, _dst_nodal); this->getDOFManager().registerDOFsDerivative("disp", 1, velocity); this->getDOFManager().registerDOFsDerivative("disp", 2, acceleration); this->getDOFManager().registerBlockedDOFs("disp", blocked); displacement.set(0.); velocity.set(0.); acceleration.set(0.); forces.set(0.); blocked.set(false); UInt global_nb_nodes = mesh.getNbGlobalNodes(); for (auto && n : arange(nb_dofs)) { auto global_id = mesh.getNodeGlobalId(n); if (global_id == (global_nb_nodes - 1)) forces(n, _x) = F; if (global_id == 0) blocked(n, _x) = true; } auto cit = this->mesh.getConnectivity(_segment_2).begin(2); auto cend = this->mesh.getConnectivity(_segment_2).end(2); auto L_it = this->initial_lengths.begin(); for (; cit != cend; ++cit, ++L_it) { const Vector & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); *L_it = std::abs(p2 - p1); } this->registerDataAccessor(*this); this->registerSynchronizer( const_cast(this->mesh.getElementSynchronizer()), SynchronizationTag::_user_1); } void assembleLumpedMass() { auto & M = this->getDOFManager().getLumpedMatrix("M"); M.clear(); this->assembleLumpedMass(_not_ghost); if (this->mesh.getNbElement(_segment_2, _ghost) > 0) this->assembleLumpedMass(_ghost); is_lumped_mass_assembled = true; } void assembleLumpedMass(const GhostType & ghost_type) { Array M(nb_dofs, 1, 0.); Array m_all_el(this->mesh.getNbElement(_segment_2, ghost_type), 2); for (auto && data : zip(make_view(this->mesh.getConnectivity(_segment_2), 2), make_view(m_all_el, 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); Real M_n = rho * A * L / 2; m_el(0) = m_el(1) = M_n; } this->getDOFManager().assembleElementalArrayLocalArray( m_all_el, M, _segment_2, ghost_type); this->getDOFManager().assembleToLumpedMatrix("disp", M, "M"); } void assembleMass() { SparseMatrix & M = this->getDOFManager().getMatrix("M"); M.clear(); 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.; } 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) { if (matrix_id == "M") { if (not is_lumped_mass_assembled) this->assembleLumpedMass(); } else { AKANTU_EXCEPTION("This solver does not know what to do with a matrix " << matrix_id); } } void assembleStiffness() { SparseMatrix & K = this->getDOFManager().getMatrix("K"); K.clear(); 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; } this->getDOFManager().assembleElementalMatricesToMatrix( "K", "disp", k_all_el, _segment_2); is_stiffness_assembled = true; } void assembleResidual() { this->getDOFManager().assembleToResidual("disp", forces); internal_forces.clear(); this->assembleResidual(_not_ghost); this->synchronize(SynchronizationTag::_user_1); this->getDOFManager().assembleToResidual("disp", internal_forces, -1.); } void assembleResidual(const GhostType & ghost_type) { Array forces_internal_el( this->mesh.getNbElement(_segment_2, ghost_type), 2); auto cit = this->mesh.getConnectivity(_segment_2, ghost_type).begin(2); auto cend = this->mesh.getConnectivity(_segment_2, ghost_type).end(2); auto f_it = forces_internal_el.begin(2); auto strain_it = this->strains.begin(); auto stress_it = this->stresses.begin(); auto L_it = this->initial_lengths.begin(); for (; cit != cend; ++cit, ++f_it, ++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; Vector & f = *f_it; f(0) = -f_n; f(1) = f_n; } this->getDOFManager().assembleElementalArrayLocalArray( forces_internal_el, internal_forces, _segment_2, ghost_type); } Real getPotentialEnergy() { Real res = 0; if (not lumped) { res = this->mulVectMatVect(this->displacement, "K", this->displacement); } else { auto strain_it = this->strains.begin(); auto stress_it = this->stresses.begin(); auto strain_end = this->strains.end(); auto L_it = this->initial_lengths.begin(); for (; strain_it != strain_end; ++strain_it, ++stress_it, ++L_it) { res += *strain_it * *stress_it * A * *L_it; } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); } return res / 2.; } Real getKineticEnergy() { Real res = 0; if (not lumped) { res = this->mulVectMatVect(this->velocity, "M", this->velocity); } else { Array & m = dynamic_cast( this->getDOFManager().getLumpedMatrix("M")); auto it = velocity.begin(); auto end = velocity.end(); auto m_it = m.begin(); for (UInt node = 0; it != end; ++it, ++m_it, ++node) { if (mesh.isLocalOrMasterNode(node)) res += *m_it * *it * *it; } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); } return res / 2.; } 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(); } Real mulVectMatVect(const Array & x, const ID & A_id, const Array & y) { Array Ay(nb_dofs); this->getDOFManager().assembleMatMulVectToArray("disp", A_id, y, Ay); Real res = 0.; for (auto && data : zip(arange(nb_dofs), make_view(Ay), make_view(x))) { res += std::get<2>(data) * std::get<1>(data) * mesh.isLocalOrMasterNode(std::get<0>(data)); } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); return res; } 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; } } } const Mesh & getMesh() const { return mesh; } UInt getSpatialDimension() const { return 1; } auto & getBlockedDOFs() { return blocked; } private: UInt nb_dofs; UInt nb_elements; bool lumped; bool is_stiffness_assembled{false}; bool is_mass_assembled{false}; bool is_lumped_mass_assembled{false}; 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; }; #endif /* __AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH__ */ } // namespace akantu diff --git a/test/test_model/test_common/test_model_solver/test_model_solver.verified b/test/test_model/test_common/test_model_solver/test_model_solver_petsc.verified similarity index 100% rename from test/test_model/test_common/test_model_solver/test_model_solver.verified rename to test/test_model/test_common/test_model_solver/test_model_solver_petsc.verified