diff --git a/packages/00_core.cmake b/packages/00_core.cmake
index bdf3e8fc3..b8e823857 100644
--- a/packages/00_core.cmake
+++ b/packages/00_core.cmake
@@ -1,432 +1,434 @@
#===============================================================================
# @file 00_core.cmake
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date creation: Mon Nov 21 2011
# @date last modification: Fri Sep 19 2014
#
# @brief package description for core
#
# @section LICENSE
#
# Copyright (©) 2010-2012, 2014 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 <http://www.gnu.org/licenses/>.
#
#===============================================================================
set(AKANTU_CORE ON CACHE INTERNAL "core package for Akantu" FORCE)
set(AKANTU_CORE_FILES
common/aka_array.cc
common/aka_array.hh
common/aka_array_tmpl.hh
common/aka_blas_lapack.hh
common/aka_circular_array.hh
common/aka_circular_array_inline_impl.cc
common/aka_common.cc
common/aka_common.hh
common/aka_common_inline_impl.cc
common/aka_csr.hh
common/aka_error.cc
common/aka_error.hh
common/aka_event_handler_manager.hh
common/aka_extern.cc
common/aka_fwd.hh
common/aka_grid_dynamic.hh
common/aka_math.cc
common/aka_math.hh
common/aka_math_tmpl.hh
common/aka_memory.cc
common/aka_memory.hh
common/aka_memory_inline_impl.cc
common/aka_random_generator.hh
common/aka_safe_enum.hh
common/aka_static_memory.cc
common/aka_static_memory.hh
common/aka_static_memory_inline_impl.cc
common/aka_static_memory_tmpl.hh
common/aka_typelist.hh
common/aka_types.hh
common/aka_visitor.hh
common/aka_voigthelper.hh
common/aka_voigthelper.cc
fe_engine/element_class.cc
fe_engine/element_class.hh
fe_engine/element_class_tmpl.hh
fe_engine/element_classes/element_class_hexahedron_8_inline_impl.cc
fe_engine/element_classes/element_class_pentahedron_6_inline_impl.cc
fe_engine/element_classes/element_class_point_1_inline_impl.cc
fe_engine/element_classes/element_class_quadrangle_4_inline_impl.cc
fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc
fe_engine/element_classes/element_class_segment_2_inline_impl.cc
fe_engine/element_classes/element_class_segment_3_inline_impl.cc
fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.cc
fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.cc
fe_engine/element_classes/element_class_triangle_3_inline_impl.cc
fe_engine/element_classes/element_class_triangle_6_inline_impl.cc
fe_engine/fe_engine.cc
fe_engine/fe_engine.hh
fe_engine/fe_engine_inline_impl.cc
fe_engine/fe_engine_template.hh
fe_engine/fe_engine_template_tmpl.hh
fe_engine/geometrical_data_tmpl.hh
fe_engine/geometrical_element.cc
fe_engine/integration_element.cc
fe_engine/integrator.hh
fe_engine/integrator_gauss.hh
fe_engine/integrator_gauss_inline_impl.cc
fe_engine/interpolation_element.cc
fe_engine/interpolation_element_tmpl.hh
fe_engine/shape_functions.hh
fe_engine/shape_functions_inline_impl.cc
fe_engine/shape_lagrange.cc
fe_engine/shape_lagrange.hh
fe_engine/shape_lagrange_inline_impl.cc
fe_engine/shape_linked.cc
fe_engine/shape_linked.hh
fe_engine/shape_linked_inline_impl.cc
fe_engine/element.hh
io/dumper/dumpable.hh
io/dumper/dumpable.cc
io/dumper/dumpable_inline_impl.hh
io/dumper/dumper_field.hh
io/dumper/dumper_material_padders.hh
io/dumper/dumper_filtered_connectivity.hh
io/dumper/dumper_element_type.hh
io/dumper/dumper_element_partition.hh
io/mesh_io.cc
io/mesh_io.hh
io/mesh_io/mesh_io_abaqus.cc
io/mesh_io/mesh_io_abaqus.hh
io/mesh_io/mesh_io_diana.cc
io/mesh_io/mesh_io_diana.hh
io/mesh_io/mesh_io_msh.cc
io/mesh_io/mesh_io_msh.hh
io/model_io.cc
io/model_io.hh
io/parser/algebraic_parser.hh
io/parser/input_file_parser.hh
io/parser/parsable.cc
io/parser/parsable.hh
io/parser/parsable_tmpl.hh
io/parser/parser.cc
io/parser/parser.hh
io/parser/parser_tmpl.hh
io/parser/cppargparse/cppargparse.hh
io/parser/cppargparse/cppargparse.cc
io/parser/cppargparse/cppargparse_tmpl.hh
mesh/element_group.cc
mesh/element_group.hh
mesh/element_group_inline_impl.cc
mesh/element_type_map.hh
mesh/element_type_map_tmpl.hh
mesh/element_type_map_filter.hh
mesh/group_manager.cc
mesh/group_manager.hh
mesh/group_manager_inline_impl.cc
mesh/mesh.cc
mesh/mesh.hh
mesh/mesh_filter.hh
mesh/mesh_data.cc
mesh/mesh_data.hh
mesh/mesh_data_tmpl.hh
mesh/mesh_inline_impl.cc
mesh/node_group.cc
mesh/node_group.hh
mesh/node_group_inline_impl.cc
mesh_utils/mesh_partition.cc
mesh_utils/mesh_partition.hh
mesh_utils/mesh_partition/mesh_partition_mesh_data.cc
mesh_utils/mesh_partition/mesh_partition_mesh_data.hh
mesh_utils/mesh_partition/mesh_partition_scotch.hh
mesh_utils/mesh_pbc.cc
mesh_utils/mesh_utils.cc
mesh_utils/mesh_utils.hh
mesh_utils/mesh_utils_inline_impl.cc
model/boundary_condition.hh
model/boundary_condition_functor.hh
model/boundary_condition_functor_inline_impl.cc
model/boundary_condition_tmpl.hh
model/integration_scheme/generalized_trapezoidal.hh
model/integration_scheme/generalized_trapezoidal_inline_impl.cc
model/integration_scheme/integration_scheme_1st_order.hh
model/integration_scheme/integration_scheme_2nd_order.hh
model/integration_scheme/newmark-beta.hh
model/integration_scheme/newmark-beta_inline_impl.cc
model/model.cc
model/model.hh
model/model_inline_impl.cc
model/solid_mechanics/material.cc
model/solid_mechanics/material.hh
model/solid_mechanics/material_inline_impl.cc
model/solid_mechanics/material_list.hh
model/solid_mechanics/material_random_internal.hh
model/solid_mechanics/material_selector.hh
model/solid_mechanics/material_selector_tmpl.hh
model/solid_mechanics/materials/internal_field.hh
model/solid_mechanics/materials/internal_field_tmpl.hh
model/solid_mechanics/materials/random_internal_field.hh
model/solid_mechanics/materials/random_internal_field_tmpl.hh
model/solid_mechanics/solid_mechanics_model.cc
model/solid_mechanics/solid_mechanics_model.hh
model/solid_mechanics/solid_mechanics_model_inline_impl.cc
model/solid_mechanics/solid_mechanics_model_mass.cc
model/solid_mechanics/solid_mechanics_model_material.cc
model/solid_mechanics/solid_mechanics_model_tmpl.hh
model/solid_mechanics/solid_mechanics_model_event_handler.hh
model/solid_mechanics/materials/plane_stress_toolbox.hh
model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh
model/solid_mechanics/materials/material_elastic.cc
model/solid_mechanics/materials/material_elastic.hh
model/solid_mechanics/materials/material_elastic_inline_impl.cc
model/solid_mechanics/materials/material_thermal.cc
model/solid_mechanics/materials/material_thermal.hh
model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc
model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh
model/solid_mechanics/materials/material_elastic_orthotropic.cc
model/solid_mechanics/materials/material_elastic_orthotropic.hh
model/solid_mechanics/materials/material_damage/material_damage.hh
model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh
model/solid_mechanics/materials/material_damage/material_marigo.cc
model/solid_mechanics/materials/material_damage/material_marigo.hh
model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc
model/solid_mechanics/materials/material_damage/material_mazars.cc
model/solid_mechanics/materials/material_damage/material_mazars.hh
model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.cc
model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc
model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh
model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.cc
model/solid_mechanics/materials/material_plastic/material_plastic.cc
model/solid_mechanics/materials/material_plastic/material_plastic.hh
model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.cc
model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc
model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh
model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.cc
model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc
model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh
solver/solver.cc
solver/solver.hh
solver/sparse_matrix.cc
solver/sparse_matrix.hh
solver/sparse_matrix_inline_impl.cc
+ solver/static_solver.hh
+ solver/static_solver.cc
synchronizer/communication_buffer.hh
synchronizer/communication_buffer_inline_impl.cc
synchronizer/data_accessor.cc
synchronizer/data_accessor.hh
synchronizer/data_accessor_inline_impl.cc
synchronizer/distributed_synchronizer.cc
synchronizer/distributed_synchronizer.hh
synchronizer/distributed_synchronizer_tmpl.hh
synchronizer/dof_synchronizer.cc
synchronizer/dof_synchronizer.hh
synchronizer/dof_synchronizer_inline_impl.cc
synchronizer/filtered_synchronizer.cc
synchronizer/filtered_synchronizer.hh
synchronizer/mpi_type_wrapper.hh
synchronizer/pbc_synchronizer.cc
synchronizer/pbc_synchronizer.hh
synchronizer/real_static_communicator.hh
synchronizer/static_communicator.cc
synchronizer/static_communicator.hh
synchronizer/static_communicator_dummy.hh
synchronizer/static_communicator_inline_impl.hh
synchronizer/synchronizer.cc
synchronizer/synchronizer.hh
synchronizer/synchronizer_registry.cc
synchronizer/synchronizer_registry.hh
)
set(AKANTU_CORE_DEB_DEPEND
libboost-dev
)
set(AKANTU_CORE_TESTS
test_csr
test_facet_element_mapping
test_facet_extraction_tetrahedron_4
test_facet_extraction_triangle_3
test_grid
test_interpolate_stress
test_local_material
test_material_damage_non_local
test_material_thermal
test_matrix
test_mesh_boundary
test_mesh_data
test_mesh_io_msh
test_mesh_io_msh_physical_names
test_mesh_partitionate_mesh_data
test_parser
test_dumper
test_pbc_tweak
test_purify_mesh
test_solid_mechanics_model_bar_traction2d
test_solid_mechanics_model_bar_traction2d_structured
test_solid_mechanics_model_bar_traction2d_structured_pbc
test_solid_mechanics_model_boundary_condition
test_solid_mechanics_model_circle_2
test_solid_mechanics_model_cube3d
test_solid_mechanics_model_cube3d_pbc
test_solid_mechanics_model_cube3d_tetra10
test_solid_mechanics_model_square
test_solid_mechanics_model_material_eigenstrain
test_static_memory
test_surface_extraction_tetrahedron_4
test_surface_extraction_triangle_3
test_vector
test_vector_iterator
test_weight
test_math
test_material_standard_linear_solid_deviatoric_relaxation
test_material_standard_linear_solid_deviatoric_relaxation_tension
test_material_plasticity
)
set(AKANTU_CORE_MANUAL_FILES
manual.sty
manual.cls
manual.tex
manual-macros.sty
manual-titlepages.tex
manual-introduction.tex
manual-gettingstarted.tex
manual-io.tex
manual-solidmechanicsmodel.tex
manual-constitutive-laws.tex
manual-lumping.tex
manual-elements.tex
manual-appendix-elements.tex
manual-appendix-materials.tex
manual-appendix-packages.tex
manual-backmatter.tex
manual-bibliography.bib
manual-bibliographystyle.bst
figures/bc_and_ic_example.pdf
figures/boundary.pdf
figures/boundary.svg
figures/dirichlet.pdf
figures/dirichlet.svg
figures/doc_wheel.pdf
figures/doc_wheel.svg
figures/dynamic_analysis.png
figures/explicit_dynamic.pdf
figures/explicit_dynamic.svg
figures/static.pdf
figures/static.svg
figures/hooke_law.pdf
figures/hot-point-1.png
figures/hot-point-2.png
figures/implicit_dynamic.pdf
figures/implicit_dynamic.svg
figures/insertion.pdf
figures/interpolate.pdf
figures/interpolate.svg
figures/problemDomain.pdf_tex
figures/problemDomain.pdf
figures/static_analysis.png
figures/stress_strain_el.pdf
figures/tangent.pdf
figures/tangent.svg
figures/vectors.pdf
figures/vectors.svg
figures/stress_strain_neo.pdf
figures/visco_elastic_law.pdf
figures/isotropic_hardening_plasticity.pdf
figures/stress_strain_visco.pdf
figures/elements/hexahedron_8.pdf
figures/elements/hexahedron_8.svg
figures/elements/quadrangle_4.pdf
figures/elements/quadrangle_4.svg
figures/elements/quadrangle_8.pdf
figures/elements/quadrangle_8.svg
figures/elements/segment_2.pdf
figures/elements/segment_2.svg
figures/elements/segment_3.pdf
figures/elements/segment_3.svg
figures/elements/tetrahedron_10.pdf
figures/elements/tetrahedron_10.svg
figures/elements/tetrahedron_4.pdf
figures/elements/tetrahedron_4.svg
figures/elements/triangle_3.pdf
figures/elements/triangle_3.svg
figures/elements/triangle_6.pdf
figures/elements/triangle_6.svg
figures/elements/xtemp.pdf
)
find_program(READLINK_COMMAND readlink)
find_program(ADDR2LINE_COMMAND addr2line)
mark_as_advanced(READLINK_COMMAND)
mark_as_advanced(ADDR2LINE_COMMAND)
include(CheckFunctionExists)
check_function_exists(clock_gettime _clock_gettime)
if(NOT _clock_gettime)
set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY ON)
else()
set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY OFF)
endif()
set(AKANTU_CORE_DOCUMENTATION
"
This package is the core engine of \\akantu. It depends on:
\\begin{itemize}
\\item A C++ compiler (\\href{http://gcc.gnu.org/}{GCC} >= 4, or \\href{https://software.intel.com/en-us/intel-compilers}{Intel}).
\\item The cross-platform, open-source \\href{http://www.cmake.org/}{CMake} build system.
\\item The \\href{http://www.boost.org/}{Boost} C++ portable libraries.
\\item The \\href{http://www.zlib.net/}{zlib} compression library.
\\end{itemize}
Under Ubuntu (14.04 LTS) the installation can be performed using the commands:
\\begin{command}
> sudo apt-get install build-essential cmake-curses-gui libboost-dev zlib1g-dev
\\end{command}
Under Mac OS X the installation requires the following steps:
\\begin{itemize}
\\item Install Xcode
\\item Install the command line tools.
\\item Install the MacPorts project which allows to automatically
download and install opensource packages.
\\end{itemize}
Then the following commands should be typed in a terminal:
\\begin{command}
> sudo port install cmake gcc48 boost
\\end{command}
")
\ No newline at end of file
diff --git a/packages/84_petsc.cmake b/packages/84_petsc.cmake
index ad0483993..dceacb9ba 100644
--- a/packages/84_petsc.cmake
+++ b/packages/84_petsc.cmake
@@ -1,41 +1,45 @@
#===============================================================================
# @file petsc.cmake
#
# @author Nicolas Richart <nicolas.richart@epfl.ch>
# @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
# @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
#
# @date Mon Nov 21 18:19:15 2011
#
# @brief package description for PETSc support
#
# @section LICENSE
#
# Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
#
#===============================================================================
add_optional_external_package(PETSc "Add PETSc support in akantu" OFF ARGS COMPONENT ARGS CXX)
add_internal_package_dependencies(petsc parallel)
set(AKANTU_PETSC_FILES
solver/petsc_matrix.hh
solver/petsc_matrix.cc
solver/petsc_matrix_inline_impl.cc
solver/solver_petsc.hh
solver/solver_petsc.cc
)
+
+set(AKANTU_PETSC_TESTS
+ test_petsc_matrix_profile
+ )
diff --git a/src/common/aka_common.cc b/src/common/aka_common.cc
index 92311f159..d5240f3aa 100644
--- a/src/common/aka_common.cc
+++ b/src/common/aka_common.cc
@@ -1,137 +1,145 @@
/**
* @file aka_common.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Jun 14 2010
* @date last modification: Mon Sep 15 2014
*
* @brief Initialization of global variables
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_static_memory.hh"
#include "static_communicator.hh"
+#include "static_solver.hh"
#include "aka_random_generator.hh"
#include "parser.hh"
#include "cppargparse.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
void initialize(int & argc, char ** & argv) {
AKANTU_DEBUG_IN();
initialize("", argc, argv);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void initialize(const std::string & input_file, int & argc, char ** & argv) {
AKANTU_DEBUG_IN();
StaticMemory::getStaticMemory();
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(argc, argv);
debug::debugger.setParallelContext(comm.whoAmI(), comm.getNbProc());
debug::initSignalHandler();
static_argparser.setParallelContext(comm.whoAmI(), comm.getNbProc());
static_argparser.setExternalExitFunction(debug::exit);
static_argparser.addArgument("--aka_input_file", "Akantu's input file",
1, cppargparse::_string, std::string());
static_argparser.addArgument("--aka_debug_level", std::string("Akantu's overall debug level") +
std::string(" (0: error, 1: exceptions, 4: warnings, 5: info, ..., 100: dump,") +
std::string(" more info on levels can be foind in aka_error.hh)"),
1,
cppargparse::_integer, int(dblWarning));
static_argparser.addArgument("--aka_print_backtrace",
"Should Akantu print a backtrace in case of error",
0,
cppargparse::_boolean, false, true);
static_argparser.parse(argc, argv, cppargparse::_remove_parsed);
std::string infile = static_argparser["aka_input_file"];
if(infile == "") infile = input_file;
debug::setDebugLevel(dblError);
if ("" != infile) {
static_parser.parse(infile);
}
long int seed;
try {
seed = static_parser.getParameter("seed", _ppsc_current_scope);
} catch (debug::Exception &) {
seed = time(NULL);
}
int dbl_level = static_argparser["aka_debug_level"];
debug::setDebugLevel(DebugLevel(dbl_level));
debug::debugger.printBacktrace(static_argparser["aka_print_backtrace"]);
seed *= (comm.whoAmI() + 1);
Rand48Generator<Real>::seed(seed);
RandGenerator<Real>::seed(seed);
AKANTU_DEBUG_INFO("Random seed set to " << seed);
+ /// initialize external solvers
+ StaticSolver::getStaticSolver().initialize(argc, argv);
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void finalize() {
AKANTU_DEBUG_IN();
+ /// finalize external solvers
+ StaticSolver::getStaticSolver().finalize();
+
if(StaticMemory::isInstantiated()) delete &(StaticMemory::getStaticMemory());
if(StaticCommunicator::isInstantiated()) {
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
comm.barrier();
delete &comm;
}
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
cppargparse::ArgumentParser & getStaticArgumentParser() {
return static_argparser;
}
/* -------------------------------------------------------------------------- */
Parser & getStaticParser() {
return static_parser;
}
/* -------------------------------------------------------------------------- */
const ParserSection & getUserParser() {
return *(static_parser.getSubSections(_st_user).first);
}
__END_AKANTU__
diff --git a/src/common/aka_config.hh.in b/src/common/aka_config.hh.in
index 44aad462b..248cf67da 100644
--- a/src/common/aka_config.hh.in
+++ b/src/common/aka_config.hh.in
@@ -1,87 +1,88 @@
/**
* @file aka_config.hh.in
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Fri Jan 13 12:34:54 2012
*
* @brief Compilation time configuration of Akantu
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_CONFIG_HH__
#define __AKANTU_AKA_CONFIG_HH__
#define AKANTU_VERSION @AKANTU_VERSION@
#define AKANTU_VERSION_MAJOR @AKANTU_MAJOR_VERSION@
#define AKANTU_VERSION_MINOR @AKANTU_MINOR_VERSION@
#define AKANTU_VERSION_PATCH @AKANTU_PATCH_VERSION@
#cmakedefine AKANTU_USE_BLAS
#cmakedefine AKANTU_USE_LAPACK
#cmakedefine AKANTU_PARALLEL
#cmakedefine AKANTU_USE_MPI
#cmakedefine AKANTU_USE_SCOTCH
#cmakedefine AKANTU_USE_PTSCOTCH
#cmakedefine AKANTU_SCOTCH_NO_EXTERN
#cmakedefine AKANTU_USE_MUMPS
+#cmakedefine AKANTU_USE_PETSC
#cmakedefine AKANTU_USE_IOHELPER
#cmakedefine AKANTU_USE_QVIEW
#cmakedefine AKANTU_USE_BLACKDYNAMITE
#cmakedefine AKANTU_USE_NLOPT
#cmakedefine AKANTU_USE_CPPARRAY
#cmakedefine AKANTU_USE_OBSOLETE_GETTIMEOFDAY
#cmakedefine AKANTU_EXTRA_MATERIALS
#cmakedefine AKANTU_DAMAGE_NON_LOCAL
#cmakedefine AKANTU_STRUCTURAL_MECHANICS
#cmakedefine AKANTU_HEAT_TRANSFER
#cmakedefine AKANTU_COHESIVE_ELEMENT
#cmakedefine AKANTU_PARALLEL_COHESIVE_ELEMENT
#cmakedefine AKANTU_IGFEM
// BOOST Section
#cmakedefine AKANTU_BOOST_CHRONO
#cmakedefine AKANTU_BOOST_SYSTEM
// Experimental part
#cmakedefine AKANTU_CORE_CXX11
// Debug tools
//#cmakedefine AKANTU_NDEBUG
#cmakedefine AKANTU_DEBUG_TOOLS
#cmakedefine READLINK_COMMAND @READLINK_COMMAND@
#cmakedefine ADDR2LINE_COMMAND @ADDR2LINE_COMMAND@
#define __aka_inline__ inline
#endif /* __AKANTU_AKA_CONFIG_HH__ */
diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc
index 52b816bf5..51e6661b8 100644
--- a/src/common/aka_error.cc
+++ b/src/common/aka_error.cc
@@ -1,362 +1,338 @@
/**
* @file aka_error.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Sep 06 2010
* @date last modification: Tue Jul 29 2014
*
* @brief handling of errors
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_config.hh"
#include "aka_error.hh"
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <csignal>
#include <execinfo.h>
#include <cxxabi.h>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <cstring>
#include <map>
#include <sys/wait.h>
#include <sys/types.h>
#include <unistd.h>
#if defined(AKANTU_USE_OBSOLETE_GETTIMEOFDAY)
# include <sys/time.h>
#else
# include <time.h>
#endif
#ifdef AKANTU_USE_MPI
-#if defined(AKANTU_CORE_CXX11)
-#if defined(__INTEL_COMPILER)
-//#pragma warning ( disable : 383 )
-#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu
-#elif (defined(__GNUC__) || defined(__GNUG__))
-# define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
-# if GCC_VERSION > 40600
-# pragma GCC diagnostic push
-# endif
-# pragma GCC diagnostic ignored "-Wliteral-suffix"
-#endif
-#endif
# include <mpi.h>
-#if defined(AKANTU_CORE_CXX11)
-#if defined(__INTEL_COMPILER)
-//#pragma warning ( disable : 383 )
-#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu
-#elif defined(__GNUG__)
-# if GCC_VERSION > 40600
-# pragma GCC diagnostic pop
-# else
-# pragma GCC diagnostic warning "-Wliteral-suffix"
-# endif
-#endif
-#endif
#endif
#define __BEGIN_AKANTU_DEBUG__ namespace akantu { namespace debug {
#define __END_AKANTU_DEBUG__ } }
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU_DEBUG__
static void printBacktraceAndExit(int sig) {
printBacktrace(sig);
debugger.exit(-50);
}
/* ------------------------------------------------------------------------ */
void initSignalHandler() {
struct sigaction action;
action.sa_handler = &printBacktraceAndExit;
sigemptyset(&(action.sa_mask));
action.sa_flags = SA_RESETHAND;
sigaction(SIGSEGV, &action, NULL);
sigaction(SIGABRT, &action, NULL);
}
/* ------------------------------------------------------------------------ */
std::string demangle(const char *symbol) {
int status;
std::string result;
char *demangled_name;
if ((demangled_name = abi::__cxa_demangle(symbol, NULL, 0, &status)) != NULL) {
result = demangled_name;
free(demangled_name);
} else {
result = symbol;
}
return result;
}
/* ------------------------------------------------------------------------ */
std::string exec(std::string cmd) {
FILE *pipe = popen(cmd.c_str(), "r");
if (!pipe) return "";
char buffer[1024];
std::string result = "";
while (!feof(pipe)) {
if (fgets(buffer, 128, pipe) != NULL)
result += buffer;
}
result = result.substr(0, result.size() - 1);
pclose(pipe);
return result;
}
/* ------------------------------------------------------------------------ */
void printBacktrace(__attribute__((unused)) int sig) {
AKANTU_DEBUG_INFO("Caught signal " << sig << "!");
#if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND)
std::string me = "";
char buf[1024];
/* The manpage says it won't null terminate. Let's zero the buffer. */
memset(buf, 0, sizeof(buf));
/* Note we use sizeof(buf)-1 since we may need an extra char for NUL. */
if (readlink("/proc/self/exe", buf, sizeof(buf) - 1))
me = std::string(buf);
std::ifstream inmaps;
inmaps.open("/proc/self/maps");
std::map <std::string, size_t> addr_map;
std::string line;
while (inmaps.good()) {
std::getline(inmaps, line);
std::stringstream sstr(line);
size_t first = line.find('-');
std::stringstream sstra(line.substr(0, first));
size_t addr; sstra >> std::hex >> addr;
std::string lib;
sstr >> lib;
sstr >> lib;
sstr >> lib;
sstr >> lib;
sstr >> lib;
sstr >> lib;
if (lib != "" && addr_map.find(lib) == addr_map.end()) {
addr_map[lib] = addr;
}
}
if (me != "") addr_map[me] = 0;
#endif
const size_t max_depth = 100;
size_t stack_depth;
void *stack_addrs[max_depth];
char **stack_strings;
size_t i;
stack_depth = backtrace(stack_addrs, max_depth);
stack_strings = backtrace_symbols(stack_addrs, stack_depth);
std::cerr << "BACKTRACE : " << stack_depth << " stack frames." << std::endl;
size_t w = size_t(std::floor(log(double(stack_depth)) / std::log(10.)) + 1);
/// -1 to remove the call to the printBacktrace function
for (i = 1; i < stack_depth; i++) {
std::cerr << std::dec << " [" << std::setw(w) << i << "] ";
std::string bt_line(stack_strings[i]);
size_t first, second;
if ((first = bt_line.find('(')) != std::string::npos && (second = bt_line.find('+')) != std::string::npos) {
std::string location = bt_line.substr(0, first);
#if defined(READLINK_COMMAND)
location = exec(std::string(BOOST_PP_STRINGIZE(READLINK_COMMAND)) + std::string(" -f ") + location);
#endif
std::string call = demangle(bt_line.substr(first + 1, second - first - 1).c_str());
size_t f = bt_line.find('[');
size_t s = bt_line.find(']');
std::string address = bt_line.substr(f + 1, s - f - 1);
std::stringstream sstra(address);
size_t addr; sstra >> std::hex >> addr;
std::cerr << location << " [" << call << "]";
#if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND)
std::map <std::string, size_t>::iterator it = addr_map.find(location);
if (it != addr_map.end()) {
std::stringstream syscom;
syscom << BOOST_PP_STRINGIZE(ADDR2LINE_COMMAND) << " 0x" << std::hex << (addr - it->second) << " -i -e " << location;
std::string line = exec(syscom.str());
std::cerr << " (" << line << ")" << std::endl;
} else {
#endif
std::cerr << " (0x" << std::hex << addr << ")" << std::endl;
#if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND)
}
#endif
} else {
std::cerr << bt_line << std::endl;
}
}
free(stack_strings);
std::cerr << "END BACKTRACE" << std::endl;
}
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
Debugger::Debugger() {
cout = &std::cerr;
level = dblWarning;
parallel_context = "";
file_open = false;
print_backtrace = false;
}
/* ------------------------------------------------------------------------ */
Debugger::~Debugger() {
if (file_open) {
dynamic_cast <std::ofstream *>(cout)->close();
delete cout;
}
}
/* ------------------------------------------------------------------------ */
void Debugger::exit(int status) {
if (status != EXIT_SUCCESS && status != -50) {
#ifndef AKANTU_NDEBUG
int * a = NULL;
*a = 1;
#endif
if(this->print_backtrace)
akantu::debug::printBacktrace(15);
}
#ifdef AKANTU_USE_MPI
if (status != EXIT_SUCCESS)
MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN);
#endif
std::exit(status); // not called when compiled with MPI due to MPI_Abort, but
// MPI_Abort does not have the noreturn attribute
}
/*------------------------------------------------------------------------- */
void Debugger::throwException(const std::string & info,
const std::string & file,
unsigned int line,
__attribute__((unused)) bool silent,
__attribute__((unused)) const std::string & location) const
throw(akantu::debug::Exception) {
#if !defined(AKANTU_NDEBUG)
if (!silent) {
printMessage("###", dblWarning, info + location);
}
#endif
debug::Exception ex(info, file, line);
throw ex;
}
/* ------------------------------------------------------------------------ */
void Debugger::printMessage(const std::string & prefix,
const DebugLevel & level,
const std::string & info) const {
if (this->level >= level) {
#if defined(AKANTU_USE_OBSOLETE_GETTIMEOFDAY)
struct timeval time;
gettimeofday(&time, NULL);
double timestamp = time.tv_sec * 1e6 + time.tv_usec; /*in us*/
#else
struct timespec time;
clock_gettime(CLOCK_REALTIME_COARSE, &time);
double timestamp = time.tv_sec * 1e6 + time.tv_nsec * 1e-3; /*in us*/
#endif
*(cout) << parallel_context
<< "{" << (unsigned int)timestamp << "} "
<< prefix << " " << info
<< std::endl;
}
}
/* ------------------------------------------------------------------------ */
void Debugger::setDebugLevel(const DebugLevel & level) {
this->level = level;
}
/* ------------------------------------------------------------------------ */
const DebugLevel & Debugger::getDebugLevel() const {
return level;
}
/* ------------------------------------------------------------------------ */
void Debugger::setLogFile(const std::string & filename) {
if (file_open) {
dynamic_cast <std::ofstream *>(cout)->close();
delete cout;
}
std::ofstream *fileout = new std::ofstream(filename.c_str());
file_open = true;
cout = fileout;
}
std::ostream & Debugger::getOutputStream() {
return *cout;
}
/* ------------------------------------------------------------------------ */
void Debugger::setParallelContext(int rank, int size) {
std::stringstream sstr;
UInt pad = std::ceil(std::log10(size));
sstr << "<" << getpid() << ">[R" << std::setfill(' ') << std::right << std::setw(pad)
<< rank << "|S" << size << "] ";
parallel_context = sstr.str();
}
void setDebugLevel(const DebugLevel & level) {
debugger.setDebugLevel(level);
}
const DebugLevel & getDebugLevel() {
return debugger.getDebugLevel();
}
/* -------------------------------------------------------------------------- */
void exit(int status) {
debugger.exit(status);
}
__END_AKANTU_DEBUG__
diff --git a/src/common/aka_extern.cc b/src/common/aka_extern.cc
index 5dd4bd8ea..b2648381c 100644
--- a/src/common/aka_extern.cc
+++ b/src/common/aka_extern.cc
@@ -1,96 +1,101 @@
/**
* @file aka_extern.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Jun 14 2010
* @date last modification: Thu Apr 03 2014
*
* @brief initialisation of all global variables
* to insure the order of creation
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_array.hh"
#include "aka_math.hh"
#include "aka_random_generator.hh"
#include "parser.hh"
#include "cppargparse.hh"
+#include "static_solver.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <limits>
/* -------------------------------------------------------------------------- */
#if defined(AKANTU_DEBUG_TOOLS)
# include "aka_debug_tools.hh"
#endif
__BEGIN_AKANTU__
/** \todo write function to get this
* values from the environment or a config file
*/
/* -------------------------------------------------------------------------- */
/* error.hpp variables */
/* -------------------------------------------------------------------------- */
namespace debug {
/// standard output for debug messages
std::ostream *_akantu_debug_cout = &std::cerr;
/// standard output for normal messages
std::ostream & _akantu_cout = std::cout;
/// parallel context used in debug messages
std::string _parallel_context = "";
Debugger debugger;
#if defined(AKANTU_DEBUG_TOOLS)
DebugElementManager element_manager;
#endif
}
/// Paser for commandline arguments
::cppargparse::ArgumentParser static_argparser;
/// Parser containing the information parsed by the input file given to initFull
Parser static_parser;
bool Parser::parser_permissive = false;
Real Math::tolerance = std::numeric_limits<Real>::epsilon();
const UInt _all_dimensions = UInt(-1);
const Array<UInt> empty_filter(0, 1, "empty_filter");
template<> long int RandGenerator<Real>::_seed = 0;
template<> long int RandGenerator<bool>::_seed = 0; // useless just defined due to a template instantiation
template<> long int RandGenerator<UInt>::_seed = 0;
template<> long int RandGenerator<Int>::_seed = 0;
template<> long int Rand48Generator<Real>::_seed = 0;
+/* -------------------------------------------------------------------------- */
+UInt StaticSolver::nb_references = 0;
+StaticSolver * StaticSolver::static_solver = NULL;
+
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.cc b/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.cc
index 22e45ab96..4aaa3aa70 100644
--- a/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.cc
+++ b/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.cc
@@ -1,222 +1,226 @@
/**
* @file material_damage_iterative.cc
* @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @date Fri Feb 14 14:54:44 2014
*
* @brief Specialization of the class material damage to damage only one gauss
* point at a time and propagate damage in a linear way. Max principal stress
* criterion is used as a failure criterion.
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_damage_iterative.hh"
#include "solid_mechanics_model.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
MaterialDamageIterative<spatial_dimension>::MaterialDamageIterative(SolidMechanicsModel & model,
- const ID & id) :
+ const ID & id) :
Material(model, id),
MaterialDamage<spatial_dimension>(model, id),
Sc("Sc", *this),
equivalent_stress("equivalent_stress", *this),
norm_max_equivalent_stress(0) {
AKANTU_DEBUG_IN();
this->registerParam("Sc", Sc, _pat_parsable, "critical stress threshold");
this->registerParam("prescribed_dam", prescribed_dam, 0.1, _pat_parsable | _pat_modifiable, "increase of damage in every step" );
this->registerParam("dam_threshold", dam_threshold, 0.8, _pat_parsable | _pat_modifiable, "damage threshold at which damage damage will be set to 1" );
this->registerParam("dam_tolerance", dam_tolerance, 0.01, _pat_parsable | _pat_modifiable, "damage tolerance to decide if quadrature point will be damageed" );
this->registerParam("max_damage", max_damage, 0.99999, _pat_parsable | _pat_modifiable, "maximum damage value" );
this->use_previous_stress = true;
this->use_previous_gradu = true;
this->Sc.initialize(1);
this->equivalent_stress.initialize(1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialDamageIterative<spatial_dimension>::computeNormalizedEquivalentStress(const Array<Real> & grad_u,
- ElementType el_type,
- GhostType ghost_type) {
+ ElementType el_type,
+ GhostType ghost_type) {
AKANTU_DEBUG_IN();
/// Vector to store eigenvalues of current stress tensor
Vector<Real> eigenvalues(spatial_dimension);
Array<Real>::const_iterator<Real> Sc_it = Sc(el_type).begin();
Array<Real>::iterator<Real> equivalent_stress_it = equivalent_stress(el_type).begin();
Array<Real>::const_matrix_iterator grad_u_it = grad_u.begin(spatial_dimension,
- spatial_dimension);
+ spatial_dimension);
Array<Real>::const_matrix_iterator grad_u_end = grad_u.end(spatial_dimension,
- spatial_dimension);
+ spatial_dimension);
Real * dam = this->damage(el_type, ghost_type).storage();
Matrix<Real> sigma(spatial_dimension, spatial_dimension);
for(;grad_u_it != grad_u_end; ++ grad_u_it) {
sigma.clear();
MaterialElastic<spatial_dimension>::computeStressOnQuad(*grad_u_it, sigma, 0.);
computeDamageAndStressOnQuad(sigma,*dam);
/// compute eigenvalues
sigma.eig(eigenvalues);
/// find max eigenvalue and normalize by tensile strength
*equivalent_stress_it = *(std::max_element(eigenvalues.storage(),
- eigenvalues.storage() + spatial_dimension)) / *(Sc_it);
+ eigenvalues.storage() + spatial_dimension)) / *(Sc_it);
++Sc_it;
++equivalent_stress_it;
++dam;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialDamageIterative<spatial_dimension>::computeAllStresses(GhostType ghost_type) {
AKANTU_DEBUG_IN();
/// reset normalized maximum equivalent stress
- if(ghost_type==_not_ghost)
+ if(ghost_type==_not_ghost)
norm_max_equivalent_stress = 0;
-
+
MaterialDamage<spatial_dimension>::computeAllStresses(ghost_type);
/// find global Gauss point with highest stress
StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator();
comm.allReduce(&norm_max_equivalent_stress, 1, _so_max);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialDamageIterative<spatial_dimension>::findMaxNormalizedEquivalentStress(ElementType el_type,
- GhostType ghost_type) {
+ GhostType ghost_type) {
AKANTU_DEBUG_IN();
if(ghost_type==_not_ghost) {
// const Array<Real> & e_stress = equivalent_stress(el_type);
// if (e_stress.begin() != e_stress.end() ) {
// Array<Real>::const_iterator<Real> equivalent_stress_it_max = std::max_element(e_stress.begin(),e_stress.end());
// /// check if max equivalent stress for this element type is greater than the current norm_max_eq_stress
// if (*equivalent_stress_it_max > norm_max_equivalent_stress)
// norm_max_equivalent_stress = *equivalent_stress_it_max;
// }
const Array<Real> & e_stress = equivalent_stress(el_type);
Array<Real>::const_iterator<Real> equivalent_stress_it = e_stress.begin();
Array<Real>::const_iterator<Real> equivalent_stress_end = e_stress.end();
Array<Real> & dam = this->damage(el_type);
Array<Real>::iterator<Real> dam_it = dam.begin();
+
for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it ) {
/// check if max equivalent stress for this element type is greater than the current norm_max_eq_stress and if the element is not already fully damaged
if (*equivalent_stress_it > norm_max_equivalent_stress && *dam_it < max_damage) {
norm_max_equivalent_stress = *equivalent_stress_it;
}
+
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialDamageIterative<spatial_dimension>::computeStress(ElementType el_type,
- GhostType ghost_type) {
+ GhostType ghost_type) {
AKANTU_DEBUG_IN();
MaterialDamage<spatial_dimension>::computeStress(el_type, ghost_type);
Real * dam = this->damage(el_type, ghost_type).storage();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
computeDamageAndStressOnQuad(sigma,*dam);
++dam;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
computeNormalizedEquivalentStress(this->gradu(el_type, ghost_type), el_type, ghost_type);
norm_max_equivalent_stress = 0;
findMaxNormalizedEquivalentStress(el_type, ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
UInt MaterialDamageIterative<spatial_dimension>::updateDamage() {
UInt nb_damaged_elements = 0;
AKANTU_DEBUG_ASSERT(prescribed_dam > 0.,
- "Your prescribed damage must be greater than zero");
+ "Your prescribed damage must be greater than zero");
if (norm_max_equivalent_stress >= 1.) {
AKANTU_DEBUG_IN();
GhostType ghost_type = _not_ghost;;
Mesh::type_iterator it = this->model->getFEEngine().getMesh().firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = this->model->getFEEngine().getMesh().lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
ElementType el_type = *it;
const Array<Real> & e_stress = equivalent_stress(el_type);
Array<Real>::const_iterator<Real> equivalent_stress_it = e_stress.begin();
Array<Real>::const_iterator<Real> equivalent_stress_end = e_stress.end();
Array<Real> & dam = this->damage(el_type);
Array<Real>::iterator<Real> dam_it = dam.begin();
for (; equivalent_stress_it != equivalent_stress_end; ++equivalent_stress_it, ++dam_it ) {
+
/// check if damage occurs
if (*equivalent_stress_it >= (1-dam_tolerance)*norm_max_equivalent_stress) {
if (*dam_it < dam_threshold)
*dam_it +=prescribed_dam;
else *dam_it = max_damage;
nb_damaged_elements += 1;
}
+
}
}
}
StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator();
comm.allReduce(&nb_damaged_elements, 1, _so_sum);
AKANTU_DEBUG_OUT();
return nb_damaged_elements;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialDamageIterative<spatial_dimension>::updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_type) {
MaterialDamage<spatial_dimension>::updateEnergies(el_type, ghost_type);
}
/* -------------------------------------------------------------------------- */
INSTANSIATE_MATERIAL(MaterialDamageIterative);
__END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.hh b/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.hh
index 6b21c4ffd..de5922931 100644
--- a/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.hh
+++ b/src/model/solid_mechanics/materials/material_damage/material_damage_iterative.hh
@@ -1,136 +1,136 @@
/**
* @file material_damage_iterative.hh
* @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @date Fri Feb 14 14:34:37 2014
*
* @brief Specialization of the class material damage to damage only one gauss
* point at a time and propagate damage in a linear way. Max principal stress
* criterion is used as a failure criterion.
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "material_damage.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__
#define __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__
__BEGIN_AKANTU__
/**
* Material damage iterative
*
* parameters in the material files :
* - Sc
*/
template<UInt spatial_dimension>
class MaterialDamageIterative : public MaterialDamage<spatial_dimension> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialDamageIterative(SolidMechanicsModel & model, const ID & id = "");
virtual ~MaterialDamageIterative() {};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// virtual void updateInternalParameters();
virtual void computeAllStresses(GhostType ghost_type = _not_ghost);
/// update internal field damage
UInt updateDamage();
/// update energies after damage has been updated
virtual void updateEnergiesAfterDamage(ElementType el_type, GhostType ghost_typ);
virtual void onBeginningSolveStep(const AnalysisMethod & method) { };
virtual void onEndSolveStep(const AnalysisMethod & method) { };
protected:
/// constitutive law for all element of a type
virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost);
///compute the equivalent stress on each Gauss point (i.e. the max prinicpal stress) and normalize it by the tensile strength
void computeNormalizedEquivalentStress(const Array<Real> & grad_u,
- ElementType el_type, GhostType ghost_type = _not_ghost);
+ ElementType el_type, GhostType ghost_type = _not_ghost);
/// find max normalized equivalent stress
void findMaxNormalizedEquivalentStress(ElementType el_type, GhostType ghost_type = _not_ghost);
inline void computeDamageAndStressOnQuad(Matrix<Real> & sigma,
- Real & dam);
+ Real & dam);
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get max normalized equivalent stress
AKANTU_GET_MACRO(NormMaxEquivalentStress, norm_max_equivalent_stress, Real);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// resistance to damage
RandomInternalField<Real> Sc;
/// internal field to store equivalent stress on each Gauss point
InternalField<Real> equivalent_stress;
/// damage increment
Real prescribed_dam;
/// maximum equivalent stress
Real norm_max_equivalent_stress;
/// deviation from max stress at which Gauss point will still get damaged
Real dam_tolerance;
/// define damage threshold at which damage will be set to 1
Real dam_threshold;
/// maximum damage value
Real max_damage;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "material_damage_iterative_inline_impl.cc"
__END_AKANTU__
#endif /* __AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH__ */
diff --git a/src/solver/petsc_matrix.cc b/src/solver/petsc_matrix.cc
index 5ad256638..e40de7b40 100644
--- a/src/solver/petsc_matrix.cc
+++ b/src/solver/petsc_matrix.cc
@@ -1,717 +1,509 @@
/**
* @file petsc_matrix.cc
* @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @date Mon Jul 21 17:40:41 2014
*
* @brief Implementation of PETSc matrix class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "petsc_matrix.hh"
+#include "static_communicator.hh"
+#include "static_communicator_mpi.hh"
+#include "mpi_type_wrapper.hh"
+#include "dof_synchronizer.hh"
/* -------------------------------------------------------------------------- */
+#include <cstring>
+#include <mpi.h>
#include <petscmat.h>
+#include <petscao.h>
+#include <petscis.h>
+#include <petscsys.h>
+
__BEGIN_AKANTU__
+struct PETScWrapper {
+ Mat mat;
+ AO ao;
+ ISLocalToGlobalMapping mapping;
+ /// pointer to the MPI communicator for PETSc commands
+ MPI_Comm communicator;
+};
+
/* -------------------------------------------------------------------------- */
-PetscMatrix::PetscMatrix(UInt size,
- const SparseMatrixType & sparse_matrix_type,
- UInt nb_degree_of_freedom,
- const ID & id,
- const MemoryID & memory_id) :
- PetscMatrix(size, sparse_matrix_type, nb_degree_of_freedom, id, memory_id) {
+PETScMatrix::PETScMatrix(UInt size,
+ const SparseMatrixType & sparse_matrix_type,
+ const ID & id,
+ const MemoryID & memory_id) :
+ SparseMatrix(size, sparse_matrix_type, id, memory_id),
+ petsc_wrapper(new PETScWrapper), d_nnz(0,1,"dnnz"), o_nnz(0,1,"onnz") {
AKANTU_DEBUG_IN();
-
- StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
- nb_proc = comm.getNbProc();
-
- UInt nb_local_nodes = model.getMesh().getNbNodes();
-
- MatCreate(&comm, mat);
-
- MatSetSizes(mat, nb_degree_of_freedom*nb_local_nodes, nb_degree_of_freedom*nb_local_nodes, size, size);
+
+ this->init();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
-PetscMatrix::PetscMatrix(const PetscMatrix & matrix,
+PETScMatrix::PETScMatrix(const SparseMatrix & matrix,
const ID & id,
const MemoryID & memory_id) :
- PetscMatrix(matrix, id, memory_id) {
+ SparseMatrix(matrix, id, memory_id),
+ petsc_wrapper(new PETScWrapper), d_nnz(0,1,"dnnz"), o_nnz(0,1,"onnz") {
AKANTU_DEBUG_IN();
- StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
- nb_proc = comm.getNbProc();
+ this->init();
- UInt nb_local_nodes = model.getMesh().getNbNodes();
- MatCreate(&comm, &mat);
- MatSetSizes(mat, nb_degree_of_freedom*nb_local_nodes, nb_degree_of_freedom*nb_local_nodes, size, size);
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+PETScMatrix::~PETScMatrix() {
+ AKANTU_DEBUG_IN();
+
+ PetscErrorCode ierr;
+ ierr = MatDestroy(&(this->petsc_wrapper->mat)); CHKERRXX(ierr);
+ delete petsc_wrapper;
AKANTU_DEBUG_OUT();
}
-
/* -------------------------------------------------------------------------- */
-PetscMatrix::~PetscMatrix() {
+void PETScMatrix::init() {
AKANTU_DEBUG_IN();
- MatDestroy(&mat);
+#if defined(AKANTU_USE_MPI)
+ StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
+ const StaticCommunicatorMPI & mpi_st_comm =
+ dynamic_cast<const StaticCommunicatorMPI &>(comm.getRealStaticCommunicator());
+ this->petsc_wrapper->communicator = mpi_st_comm.getMPITypeWrapper().getMPICommunicator();
+#else
+ this->petsc_wrapper->communicator = PETSC_COMM_SELF;
+#endif
+
+ PetscErrorCode ierr;
+
+ ierr = MatCreate(this->petsc_wrapper->communicator, &(this->petsc_wrapper->mat)); CHKERRXX(ierr);
+ /// set matrix type
+ ierr = MatSetType(this->petsc_wrapper->mat, MATAIJ); CHKERRXX(ierr);
+
+
+ if (this->sparse_matrix_type==_symmetric)
+ /// set flag for symmetry to enable ICC/Cholesky preconditioner
+ ierr = MatSetOption(this->petsc_wrapper->mat, MAT_SYMMETRIC, PETSC_TRUE); CHKERRXX(ierr);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
-PetscMatrix::createAkantuToPetscMap(Mesh & mesh) {
+void PETScMatrix::resize(const DOFSynchronizer & dof_synchronizer) {
AKANTU_DEBUG_IN();
- StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
- Int nb_proc = comm.getNbProc();
- Int rank = comm.whoAmI();
+ PetscErrorCode ierr;
this->dof_synchronizer = &const_cast<DOFSynchronizer &>(dof_synchronizer);
- //initialize vector to store the number of local and master nodes for each processor
- Vector<UInt> master_local_ndofs_per_proc(nb_proc);
+ /// find the number of dofs corresponding to master or local nodes
+ UInt nb_dofs = dof_synchronizer.getNbDOFs();
+ UInt nb_local_master_dofs = 0;
- //find the number of dofs corresponding to master or local nodes and store their akantu global dof number
- Real nb_dofs = dof_sychronizer.getNbDOFs();
- Real nb_local_master_dofs = 0;
+ /// create array to store the global equation number of all local and master dofs
+ Array<Int> local_master_eq_nbs(nb_dofs);
+ Array<Int>::scalar_iterator it_eq_nb = local_master_eq_nbs.begin();
- Array<Real> local_master_dofs(nb_dofs);
- Array<Real>::scalar_iterator it_dof = local_master_dofs.begin();
+ Array<Int> dof_types = dof_synchronizer.getDOFTypes();
+ Array<Int>::scalar_iterator it_dof_type = dof_types.begin();
- Array<Real> dof_types = dof_sychronizer.getDOFTypes();
- Array<Real>::scalar_iterator it_dof_type = dof_types.begin();
+ /// get the pointer to the gloabl equation number array
+ Int * eq_nb_val = dof_synchronizer.getGlobalDOFEquationNumbers().storage();
- for (UInt i = 0; i <nb_dofs; ++i; ++it_dof_type) {
+ for (UInt i = 0; i <nb_dofs; ++i, ++it_dof_type) {
if (*it_dof_type == -2 || *it_dof_type == -1 ) {
- *it_dof = dof_sychronizer.dof_global_ids(i);
- ++it_dof;
+ /// get for the local dof the global equation number
+ *it_eq_nb = eq_nb_val[i];
+ ++it_eq_nb;
++nb_local_master_dofs;
}
}
- local_master_dofs.resize(nb_local_master_dofs);
-
- // store the nb of master and local dofs on each processor
- master_local_ndofs_per_proc(rank) = nb_local_master_dofs;
-
- //exchange the information among all processors
- comm.allGather(local_master_ndofs_per_proc.storage(), 1);
-
- //each processor creates a map for his akantu global dofs to the corresponding petsc global dofs
- UInt petsc_start_index = 0;
+ local_master_eq_nbs.resize(nb_local_master_dofs);
- for (UInt i = 0; i < rank; ++i) {
- petsc_start_index += local_master_ndofs_per_proc(rank);
- }
+ /// set the local size
+ this->local_size = nb_local_master_dofs;
- Int * local_master_dofs_ptr = local_master_dofs.storage();
- AOCreateBasic(&comm, nb_local_master_dofs,index_akantu, index_petsc, &ao);
+ /// resize PETSc matrix
+#if defined(AKANTU_USE_MPI)
+ ierr = MatSetSizes(this->petsc_wrapper->mat, this->size, this->size, this->local_size, this->local_size); CHKERRXX(ierr);
+#else
+ ierr = MatSetSizes(this->petsc_wrapper->mat, this->local_size, this->local_size); CHKERRXX(ierr);
+#endif
+
+ /// create mapping from akantu global numbering to petsc global numbering
+ this->createGlobalAkantuToPETScMap(local_master_eq_nbs.storage());
+
AKANTU_DEBUG_OUT();
}
+
/* -------------------------------------------------------------------------- */
-void PetscMatrix::buildProfile(const Mesh & mesh, const DOFSynchronizer & dof_synchronizer) {
+void PETScMatrix::createGlobalAkantuToPETScMap(Int* local_master_eq_nbs_ptr) {
AKANTU_DEBUG_IN();
- // if(irn_jcn_to_k) delete irn_jcn_to_k;
- // irn_jcn_to_k = new std::map<std::pair<UInt, UInt>, UInt>;
- clearProfile();
-
- this->dof_synchronizer = &const_cast<DOFSynchronizer &>(dof_synchronizer);
+ PetscErrorCode ierr;
+
+ StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
+ UInt rank = comm.whoAmI();
- coordinate_list_map::iterator irn_jcn_k_it;
+ //initialize vector to store the number of local and master nodes that are assigned to each processor
+ Vector<UInt> master_local_ndofs_per_proc(nb_proc);
- Int * eq_nb_val = dof_synchronizer.getGlobalDOFEquationNumbers().storage();
+ /// store the nb of master and local dofs on each processor
+ master_local_ndofs_per_proc(rank) = this->local_size;
+
- Mesh::type_iterator it = mesh.firstType(mesh.getSpatialDimension(), _not_ghost, _ek_not_defined);
- Mesh::type_iterator end = mesh.lastType (mesh.getSpatialDimension(), _not_ghost, _ek_not_defined);
- for(; it != end; ++it) {
- UInt nb_element = mesh.getNbElement(*it);
- UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
- UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom;
- UInt * conn_val = mesh.getConnectivity(*it, _not_ghost).storage();
- Int * local_eq_nb_val = new Int[nb_degree_of_freedom * nb_nodes_per_element];
+ /// exchange the information among all processors
+ comm.allGather(master_local_ndofs_per_proc.storage(), 1);
+ /// each processor creates a map for his akantu global dofs to the corresponding petsc global dofs
- for (UInt e = 0; e < nb_element; ++e) {
- Int * tmp_local_eq_nb_val = local_eq_nb_val;
- for (UInt i = 0; i < nb_nodes_per_element; ++i) {
- UInt n = conn_val[i];
- for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
- *tmp_local_eq_nb_val++ = eq_nb_val[n * nb_degree_of_freedom + d];
- }
- // memcpy(tmp_local_eq_nb_val, eq_nb_val + n * nb_degree_of_freedom, nb_degree_of_freedom * sizeof(Int));
- // tmp_local_eq_nb_val += nb_degree_of_freedom;
- }
-
- for (UInt i = 0; i < size_mat; ++i) {
- UInt c_irn = local_eq_nb_val[i];
- if(c_irn < size) {
- UInt j_start = (sparse_matrix_type == _symmetric) ? i : 0;
- for (UInt j = j_start; j < size_mat; ++j) {
- UInt c_jcn = local_eq_nb_val[j];
- if(c_jcn < size) {
- MatSetValue(this->mat,1,c_irn, 1, c_jcn, INSERT_VALUES);
-
- if (irn_jcn_k_it == irn_jcn_k.end()) {
- irn_jcn_k[irn_jcn] = nb_non_zero;
- irn.push_back(c_irn + 1);
- jcn.push_back(c_jcn + 1);
- nb_non_zero++;
- }
- }
- }
- }
- }
- conn_val += nb_nodes_per_element;
- }
+ /// determine the PETSc-index for the first dof on each processor
+ UInt petsc_start_index = 0;
- delete [] local_eq_nb_val;
+ for (UInt i = 0; i < rank; ++i) {
+ petsc_start_index += master_local_ndofs_per_proc(rank);
}
- /// for pbc @todo correct it for parallel
- if(StaticCommunicator::getStaticCommunicator().getNbProc() == 1) {
- for (UInt i = 0; i < size; ++i) {
- KeyCOO irn_jcn = key(i, i);
- irn_jcn_k_it = irn_jcn_k.find(irn_jcn);
- if(irn_jcn_k_it == irn_jcn_k.end()) {
- irn_jcn_k[irn_jcn] = nb_non_zero;
- irn.push_back(i + 1);
- jcn.push_back(i + 1);
- nb_non_zero++;
- }
- }
+ /// create array for petsc ordering
+ Array<Int> petsc_dofs(this->local_size);
+ Array<Int>::scalar_iterator it_petsc_dofs = petsc_dofs.begin();
+
+ for (UInt i = petsc_start_index; i < petsc_start_index + this->local_size; ++i, ++it_petsc_dofs) {
+ *it_petsc_dofs = i;
}
- a.resize(nb_non_zero);
+ ierr = AOCreateBasic(this->petsc_wrapper->communicator, this->local_size, local_master_eq_nbs_ptr, petsc_dofs.storage(), &(this->petsc_wrapper->ao)); CHKERRXX(ierr);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
-void PetscMatrix::applyBoundaryNormal(Array<bool> & boundary_normal, Array<Real> & EulerAngles, Array<Real> & rhs, const Array<Real> & matrix, Array<Real> & rhs_rotated) {
- AKANTU_DEBUG_IN();
-
-
- const UInt dim = nb_degree_of_freedom;
- const UInt nb_nodes = boundary_normal.getSize();
- Matrix<Real> Ti(dim, dim);
- Matrix<Real> Tj(dim, dim);
- Matrix<Real> small_matrix(dim, dim);
- Matrix<Real> Ti_small_matrix(dim, dim);
- Matrix<Real> Ti_small_matrix_Tj(dim, dim);
- Matrix<Real> small_rhs(dim, 1);
- Matrix<Real> Ti_small_rhs(dim, 1);
- //Transformation matrix based on euler angles X_1Y_2Z_3
-
- const DOFSynchronizer::GlobalEquationNumberMap & local_eq_num_to_global = dof_synchronizer->getGlobalEquationNumberToLocal();
-
- Int * irn_val = irn.storage();
- Int * jcn_val = jcn.storage();
- Int * eq_nb_val = dof_synchronizer->getGlobalDOFEquationNumbers().storage();
- Int * eq_nb_rhs_val = dof_synchronizer->getLocalDOFEquationNumbers().storage();
-
- Array<bool> * nodes_rotated = new Array<bool > (nb_nodes, dim, "nodes_rotated");
- nodes_rotated->clear();
-
- Array<Int> * local_eq_nb_val_node_i = new Array<Int > (1, nb_degree_of_freedom, "local_eq_nb_val_node_i");
- local_eq_nb_val_node_i->clear();
- Array<Int> * local_eq_nb_val_node_j = new Array<Int > (1, nb_degree_of_freedom, "local_eq_nb_val_node_j");
- local_eq_nb_val_node_j->clear();
-
- for (UInt i = 0; i < nb_non_zero; ++i) {
- UInt ni = local_eq_num_to_global.find(*irn_val - 1)->second;
- UInt node_i = (ni - ni % nb_degree_of_freedom) / nb_degree_of_freedom;
- UInt nj = local_eq_num_to_global.find(*jcn_val - 1)->second;
- UInt node_j = (nj - nj % nb_degree_of_freedom) / nb_degree_of_freedom;
- bool constrain_ij = false;
- for (UInt j = 0; j < dim; j++){
- if ( (boundary_normal(node_i, j) || boundary_normal(node_j, j)) ) {
- constrain_ij = true;
- break;
- }
- }
-
- if(constrain_ij){
- if(dim==2){
- Real Theta=EulerAngles(node_i,0);
- Ti(0,0)=cos(Theta);
- Ti(0,1)=-sin(Theta);
- Ti(1,1)=cos(Theta);
- Ti(1,0)=sin(Theta);
-
- Theta=EulerAngles(node_j,0);
- Tj(0,0)=cos(Theta);
- Tj(0,1)=-sin(Theta);
- Tj(1,1)=cos(Theta);
- Tj(1,0)=sin(Theta);
- }
- else if(dim==3){
- Real Theta_x=EulerAngles(node_i,0);
- Real Theta_y=EulerAngles(node_i,1);
- Real Theta_z=EulerAngles(node_i,2);
-
- Ti(0,0)=cos(Theta_y)*cos(Theta_z);
- Ti(0,1)=-cos(Theta_y)*sin(Theta_z);
- Ti(0,2)=sin(Theta_y);
- Ti(1,0)=cos(Theta_x)*sin(Theta_z)+cos(Theta_z)*sin(Theta_x)*sin(Theta_y);
- Ti(1,1)=cos(Theta_x)*cos(Theta_z)-sin(Theta_x)*sin(Theta_y)*sin(Theta_z);
- Ti(1,2)=-cos(Theta_y)*sin(Theta_x);
- Ti(2,0)=sin(Theta_x)*sin(Theta_z)-cos(Theta_x)*cos(Theta_z)*sin(Theta_y);
- Ti(2,1)=cos(Theta_z)*sin(Theta_x)+cos(Theta_x)*sin(Theta_y)*sin(Theta_z);
- Ti(2,2)=cos(Theta_x)*cos(Theta_y);
-
- Theta_x=EulerAngles(node_j,0);
- Theta_y=EulerAngles(node_j,1);
- Theta_z=EulerAngles(node_j,2);
-
- Tj(0,0)=cos(Theta_y)*cos(Theta_z);
- Tj(0,1)=-cos(Theta_y)*sin(Theta_z);
- Tj(0,2)=sin(Theta_y);
- Tj(1,0)=cos(Theta_x)*sin(Theta_z)+cos(Theta_z)*sin(Theta_x)*sin(Theta_y);
- Tj(1,1)=cos(Theta_x)*cos(Theta_z)-sin(Theta_x)*sin(Theta_y)*sin(Theta_z);
- Tj(1,2)=-cos(Theta_y)*sin(Theta_x);
- Tj(2,0)=sin(Theta_x)*sin(Theta_z)-cos(Theta_x)*cos(Theta_z)*sin(Theta_y);
- Tj(2,1)=cos(Theta_z)*sin(Theta_x)+cos(Theta_x)*sin(Theta_y)*sin(Theta_z);
- Tj(2,2)=cos(Theta_x)*cos(Theta_y);
- }
- for (UInt j = 0; j < nb_degree_of_freedom; ++j){
- local_eq_nb_val_node_i->storage()[j] = eq_nb_val[node_i * nb_degree_of_freedom + j];
- local_eq_nb_val_node_j->storage()[j] = eq_nb_val[node_j * nb_degree_of_freedom + j];
- }
- for (UInt j = 0; j < nb_degree_of_freedom; ++j) {
- for (UInt k = 0; k < nb_degree_of_freedom; ++k) {
- KeyCOO jcn_irn = key(local_eq_nb_val_node_i->storage()[j], local_eq_nb_val_node_j->storage()[k]);
- coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn);
- small_matrix(j, k) = matrix.storage()[irn_jcn_k_it->second];
- }
- small_rhs(j,0)=rhs.storage()[eq_nb_rhs_val[node_i*dim+j]];
- }
-
- Ti_small_matrix.mul < false, false > (Ti, small_matrix);
- Ti_small_matrix_Tj.mul < false, true> (Ti_small_matrix, Tj);
- Ti_small_rhs.mul<false, false>(Ti, small_rhs);
-
- /*for (UInt j = 0; j < nb_degree_of_freedom; ++j) {
- for (UInt k = 0; k < nb_degree_of_freedom; ++k) {
- KeyCOO jcn_irn = key(local_eq_nb_val_node_i->storage()[j], local_eq_nb_val_node_j->storage()[k]);
- coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn);
- a.storage()[irn_jcn_k_it->second] = T_small_matrix_T(j,k);
- }
- if(node_constrain==node_i && !( nodes_rotated->storage()[node_i]))
- rhs.storage()[eq_nb_rhs_val[node_i*dim+j]]=T_small_rhs(j,0);
- }*/
- KeyCOO jcn_irn = key(ni, nj);
- coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn);
- a.storage()[irn_jcn_k_it->second] = Ti_small_matrix_Tj(ni % nb_degree_of_freedom, nj % nb_degree_of_freedom);
- //this->saveMatrix("stiffness_partial.out");
- if (!(nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]])) {
- rhs_rotated.storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = Ti_small_rhs(ni % nb_degree_of_freedom, 0);
- nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = true;
- }
- }
- else{
- if (!(nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]])) {
- rhs_rotated.storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = rhs.storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]];
- nodes_rotated->storage()[eq_nb_rhs_val[node_i * dim + ni % nb_degree_of_freedom]] = true;
- }
- }
- irn_val++;
- jcn_val++;
- }
+// void PETScMatrix::createLocalAkantuToPETScMap(const DOFSynchronizer & dof_synchronizer) {
+// AKANTU_DEBUG_IN();
- delete local_eq_nb_val_node_i;
- delete local_eq_nb_val_node_j;
- delete nodes_rotated;
- this->saveMatrix("stiffness_rotated.out");
- applyBoundary(boundary_normal);
-
- /*Real norm = 0;
- UInt n_zeros = 0;
- for (UInt j = 0; j < nb_degree_of_freedom; ++j) {
- norm += Normal.storage()[j] * Normal.storage()[j];
- if (std::abs(Normal.storage()[j]) < Math::getTolerance())
- n_zeros++;
- }
- norm = sqrt(norm);
- AKANTU_DEBUG_ASSERT(norm > Math::getTolerance(), "The normal is not right");
- for (UInt j = 0; j < nb_degree_of_freedom; ++j)
- Normal.storage()[j] /= norm;
-
- if (n_zeros == nb_degree_of_freedom - 1) {
- UInt nb_nodes = boundary_normal.getSize();
- for (UInt i = 0; i < nb_nodes; ++i) {
- if (boundary_normal(i, 0))
- for (UInt j = 0; j < nb_degree_of_freedom; ++j)
- boundary(i, j) += Normal(0, j);
- }
- } else {
-
- const DOFSynchronizer::GlobalEquationNumberMap & local_eq_num_to_global = dof_synchronizer->getGlobalEquationNumberToLocal();
-
- Int * irn_val = irn.storage();
- Int * eq_nb_val = dof_synchronizer->getGlobalDOFEquationNumbers().storage();
- Array<Int> * local_eq_nb_val = new Array<Int > (1, nb_degree_of_freedom, "local_eq_nb_val");
- local_eq_nb_val->clear();
- UInt nb_nodes = boundary_normal.getSize();
- Array<bool> * node_set = new Array<bool > (nb_nodes, 1, "node_set");
- node_set->clear();
-
- for (UInt i = 0; i < nb_non_zero; ++i) {
- UInt ni = local_eq_num_to_global.find(*irn_val - 1)->second;
- UInt node_i = (ni - ni % nb_degree_of_freedom) / nb_degree_of_freedom;
- if (boundary_normal.storage()[ni]) {
- if ((!number.storage()[ni]) && (!node_set->storage()[node_i])) {
- UInt normal_component_i = ni % nb_degree_of_freedom; //DOF of node node_i to be removed
- if (std::abs(Normal.storage()[normal_component_i]) > Math::getTolerance()) {
- for (UInt j = 0; j < nb_degree_of_freedom; ++j)
- local_eq_nb_val->storage()[j] = eq_nb_val[node_i * nb_degree_of_freedom + j];
-
- UInt DOF_remove = local_eq_nb_val->storage()[normal_component_i];
- KeyCOO jcn_irn = key(DOF_remove, DOF_remove);
- coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn);
-
- Real aii = a.storage()[irn_jcn_k_it->second];
- Real normal_constant = (1 - aii) / (Normal.storage()[normal_component_i] * Normal.storage()[normal_component_i]);
-
- for (UInt j = 0; j < nb_degree_of_freedom; ++j) {
- UInt line_j = local_eq_nb_val->storage()[j];
- for (UInt k = 0; k < nb_degree_of_freedom; ++k) {
- UInt column_k = local_eq_nb_val->storage()[k];
- jcn_irn = key(line_j, column_k);
- if ((line_j == column_k) && (line_j == DOF_remove))
- a.storage()[irn_jcn_k_it->second] = 1;
- else
- a.storage()[irn_jcn_k_it->second] += Normal.storage()[j] * Normal.storage()[k] * normal_constant;
- }
- }
-
- number.storage()[ni]++;
- node_set->storage()[node_i] = false;
- }
- }
- }
- irn_val++;
- }
-
- delete local_eq_nb_val;
- delete node_set;
- }*/
- AKANTU_DEBUG_OUT();
-}
+// AKANTU_DEBUG_ASSERT(this->petsc_wrapper->ao != NULL,
+// "You should first create a mapping from the global"
+// << " Akantu numbering to the global PETSc numbering");
-/* -------------------------------------------------------------------------- */
-void PetscMatrix::applyBoundary(const Array<bool> & boundary, Real block_val) {
- AKANTU_DEBUG_IN();
+// this->dof_synchronizer = &const_cast<DOFSynchronizer &>(dof_synchronizer);
- const DOFSynchronizer::GlobalEquationNumberMap & local_eq_num_to_global = dof_synchronizer->getGlobalEquationNumberToLocal();
- Int * irn_val = irn.storage();
- Int * jcn_val = jcn.storage();
- Real * a_val = a.storage();
-
- for (UInt i = 0; i < nb_non_zero; ++i) {
- UInt ni = local_eq_num_to_global.find(*irn_val - 1)->second;
- UInt nj = local_eq_num_to_global.find(*jcn_val - 1)->second;
- if(boundary.storage()[ni] || boundary.storage()[nj]) {
- if (*irn_val != *jcn_val) {
- *a_val = 0;
- } else {
- if(dof_synchronizer->getDOFTypes()(ni) >= 0) {
- *a_val = 0;
- } else {
- *a_val = block_val;
- }
- }
- }
- irn_val++; jcn_val++; a_val++;
- }
+// /// get the number of dofs
+// Int nb_dofs = dof_synchronizer.getNbDOFs();
- AKANTU_DEBUG_OUT();
-}
+// /// get the global equation numbers for each node
+// Array<Int> global_dof_equation_numbers = dof_synchronizer.getGlobalDOFEquationNumbers();
+
+// /// map the global dof equation numbers to the corresponding PETSc ordering
+// AOApplicationToPETSc(this->petsc_wrapper->ao, nb_dofs,
+// global_dof_equation_numbers.storage());
+
+// /// create the mapping from the local Akantu ordering to the global PETSc ordering
+// ISLocalToGlobalMappingCreate(this->petsc_wrapper->communicator,
+// 1, nb_dofs, global_dof_equation_numbers.storage(),
+// PETSC_COPY_VALUES, mapping);
+
+// AKANTU_DEBUG_OUT();
+// }
/* -------------------------------------------------------------------------- */
-void PetscMatrix::removeBoundary(const Array<bool> & boundary) {
+void PETScMatrix::buildProfile(const Mesh & mesh, const DOFSynchronizer & dof_synchronizer, UInt nb_degree_of_freedom) {
AKANTU_DEBUG_IN();
- if(irn_save) delete irn_save;
- if(jcn_save) delete jcn_save;
+ //clearProfile();
- irn_save = new Array<Int>(irn, true);
- jcn_save = new Array<Int>(jcn, true);
+ PetscErrorCode ierr;
- UInt n = boundary.getSize()*boundary.getNbComponent();
+ /// resize arrays to store the number of nonzeros in each row
+ this->d_nnz.resize(this->local_size);
+ this->o_nnz.resize(this->local_size);
- UInt * perm = new UInt[n];
+ /// set arrays to zero everywhere
+ this->d_nnz.set(0);
+ this->o_nnz.set(0);
- size_save = size;
- size = 0;
- for (UInt i = 0; i < n; ++i) {
- if(!boundary.storage()[i]) {
- perm[i] = size;
- // std::cout << "perm["<< i <<"] = " << size << std::endl;
- size++;
- }
- }
+ this->dof_synchronizer = &const_cast<DOFSynchronizer &>(dof_synchronizer);
- for (UInt i = 0; i < nb_non_zero;) {
- if(boundary.storage()[irn(i) - 1] || boundary.storage()[jcn(i) - 1]) {
- irn.erase(i);
- jcn.erase(i);
- a.erase(i);
- nb_non_zero--;
- } else {
- irn(i) = perm[irn(i) - 1] + 1;
- jcn(i) = perm[jcn(i) - 1] + 1;
- i++;
- }
- }
- delete [] perm;
+ /// get the types for each local dof
+ const Array<Int> & dof_types = dof_synchronizer.getDOFTypes();
+ Int* dof_type_ptr = dof_types.storage();
+ // if(irn_jcn_to_k) delete irn_jcn_to_k;
+ // irn_jcn_to_k = new std::map<std::pair<UInt, UInt>, UInt>;
+
- AKANTU_DEBUG_OUT();
-}
+ this->dof_synchronizer = &const_cast<DOFSynchronizer &>(dof_synchronizer);
-/* -------------------------------------------------------------------------- */
-void PetscMatrix::restoreProfile() {
- AKANTU_DEBUG_IN();
+ coordinate_list_map::iterator irn_jcn_k_it;
- irn.resize(irn_save->getSize());
- jcn.resize(jcn_save->getSize());
+ Int * eq_nb_val = dof_synchronizer.getGlobalDOFEquationNumbers().storage();
- nb_non_zero = irn.getSize();
- a.resize(nb_non_zero);
- size = size_save;
+ /// Loop over all the ghost types
+ for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) {
+ const GhostType & ghost_type = *gt;
+ /// Loop over all element types
+ Mesh::type_iterator it = mesh.firstType(mesh.getSpatialDimension(), ghost_type, _ek_not_defined);
+ Mesh::type_iterator end = mesh.lastType (mesh.getSpatialDimension(), ghost_type, _ek_not_defined);
- memcpy(irn.storage(), irn_save->storage(), irn.getSize()*sizeof(Int));
- memcpy(jcn.storage(), jcn_save->storage(), jcn.getSize()*sizeof(Int));
+ for(; it != end; ++it) {
+ UInt nb_element = mesh.getNbElement(*it);
+ UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
+ UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom;
- delete irn_save; irn_save = NULL;
- delete jcn_save; jcn_save = NULL;
+ UInt * conn_val = mesh.getConnectivity(*it, _not_ghost).storage();
+ Int * global_eq_nb_val = new Int[nb_degree_of_freedom * nb_nodes_per_element];
+ Int * local_eq_nb_val = new Int[nb_degree_of_freedom * nb_nodes_per_element];
- AKANTU_DEBUG_OUT();
-}
-/* -------------------------------------------------------------------------- */
-void PetscMatrix::saveProfile(const std::string & filename) const {
- AKANTU_DEBUG_IN();
+ for (UInt e = 0; e < nb_element; ++e) {
+ Int * tmp_global_eq_nb_val = global_eq_nb_val;
+ Int * tmp_local_eq_nb_val = local_eq_nb_val;
- std::ofstream outfile;
- outfile.open(filename.c_str());
+ for (UInt i = 0; i < nb_nodes_per_element; ++i) {
+ UInt n = conn_val[i];
+ for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
+ *tmp_global_eq_nb_val++ = eq_nb_val[n * nb_degree_of_freedom + d];
+ *tmp_local_eq_nb_val++ = n * nb_degree_of_freedom + d;
+ }
+ // memcpy(tmp_global_eq_nb_val, eq_nb_val + n * nb_degree_of_freedom, nb_degree_of_freedom * sizeof(Int));
+ // tmp_global_eq_nb_val += nb_degree_of_freedom;
+ }
- outfile << "%%MatrixMarket matrix coordinate pattern";
+ for (UInt i = 0; i < size_mat; ++i) {
+ if (dof_type_ptr[local_eq_nb_val[i]] == -2 || dof_type_ptr[local_eq_nb_val[i]] == -1 ) {
+ UInt c_irn = global_eq_nb_val[i];
+ if(c_irn < size ) {
+ //UInt j_start = (sparse_matrix_type == _symmetric) ? i : 0;
+ for (UInt j = 0; j < size_mat; ++j) {
+ UInt c_jcn = global_eq_nb_val[j];
+ if(c_jcn < size) {
+ KeyCOO irn_jcn = key(c_irn, c_jcn);
+ irn_jcn_k_it = irn_jcn_k.find(irn_jcn);
+
+ if (irn_jcn_k_it == irn_jcn_k.end()) {
+ if (dof_type_ptr[local_eq_nb_val[j]] == -2 || dof_type_ptr[local_eq_nb_val[j]] == -1 )
+ this->d_nnz(local_eq_nb_val[i]) += 1;
+ else
+ this->o_nnz(local_eq_nb_val[i]) += 1;
+ }
+ }
+ }
+ }
+ }
+ }
+ conn_val += nb_nodes_per_element;
+ }
- if(sparse_matrix_type == _symmetric) outfile << " symmetric";
- else outfile << " general";
- outfile << std::endl;
+ delete [] global_eq_nb_val;
+ delete [] local_eq_nb_val;
+ }
- UInt m = size;
- outfile << m << " " << m << " " << nb_non_zero << std::endl;
+ // /// for pbc @todo correct it for parallel
+ // if(StaticCommunicator::getStaticCommunicator().getNbProc() == 1) {
+ // for (UInt i = 0; i < size; ++i) {
+ // KeyCOO irn_jcn = key(i, i);
+ // irn_jcn_k_it = irn_jcn_k.find(irn_jcn);
+ // if(irn_jcn_k_it == irn_jcn_k.end()) {
+ // irn_jcn_k[irn_jcn] = nb_non_zero;
+ // irn.push_back(i + 1);
+ // jcn.push_back(i + 1);
+ // nb_non_zero++;
+ // }
+ // }
+ // }
+ }
- for (UInt i = 0; i < nb_non_zero; ++i) {
- outfile << irn.storage()[i] << " " << jcn.storage()[i] << " 1" << std::endl;
+ /// map arrays of nonzeros per row from akantu global to petsc global numbering
+ ierr = AOApplicationToPetsc(this->petsc_wrapper->ao,
+ local_size, d_nnz.storage()); CHKERRXX(ierr);
+ ierr = AOApplicationToPetsc(this->petsc_wrapper->ao,
+ local_size, o_nnz.storage()); CHKERRXX(ierr);
+
+ // std::string mat_type;
+ // mat_type.resize(20, 'a');
+ // std::cout << "MatType: " << mat_type << std::endl;
+ // const char * mat_type_ptr = mat_type.c_str();
+
+ MatType type;
+ MatGetType(this->petsc_wrapper->mat, &type);
+
+ // PetscTypeCompare((PetscObject)(this->petsc_wrapper->mat),MATSEQAIJ,&sametype);
+ //ierr = PetscTypeCompare(pet, MATSEQAIJ, &sametype); CHKERRXX(ierr);
+
+ /// build profile:
+ if (strcmp(type, MATSEQAIJ) == 0) {
+ ierr = MatSeqAIJSetPreallocation(this->petsc_wrapper->mat,
+ 0, d_nnz.storage()); CHKERRXX(ierr);
+ } else if ((strcmp(type, MATMPIAIJ) == 0)) {
+ ierr = MatMPIAIJSetPreallocation(this->petsc_wrapper->mat,
+ 0, d_nnz.storage(), 0,
+ o_nnz.storage()); CHKERRXX(ierr);
+ } else {
+ AKANTU_DEBUG_ERROR("The type " << type << " of PETSc matrix is not handled by"
+ << " akantu in the preallocation step");
}
- outfile.close();
+ /// assemble matrix
+ MatAssemblyBegin(this->petsc_wrapper->mat, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(this->petsc_wrapper->mat, MAT_FINAL_ASSEMBLY);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
-void PetscMatrix::saveMatrix(const std::string & filename) const {
+void PETScMatrix::saveMatrix(const std::string & filename) const {
AKANTU_DEBUG_IN();
- std::ofstream outfile;
- outfile.precision(std::numeric_limits<Real>::digits10);
-
- outfile.open(filename.c_str());
-
- outfile << "%%MatrixMarket matrix coordinate real";
-
- if(sparse_matrix_type == _symmetric) outfile << " symmetric";
- else outfile << " general";
- outfile << std::endl;
-
- outfile << size << " " << size << " " << nb_non_zero << std::endl;
-
- for (UInt i = 0; i < nb_non_zero; ++i) {
- outfile << irn(i) << " " << jcn(i) << " " << a(i) << std::endl;
- }
+ PetscErrorCode ierr;
+
+ /// create Petsc viewer
+ PetscViewer viewer;
+ ierr = PetscViewerASCIIOpen(this->petsc_wrapper->communicator, filename.c_str(), &viewer); CHKERRXX(ierr);
+ /// save the matrix
+ ierr = MatView(this->petsc_wrapper->mat, viewer); CHKERRXX(ierr);
- outfile.close();
+ /// destroy the viewer
+ ierr = PetscViewerDestroy(&viewer); CHKERRXX(ierr);
AKANTU_DEBUG_OUT();
}
-/* -------------------------------------------------------------------------- */
-Array<Real> & operator*=(Array<Real> & vect, const PetscMatrix & mat) {
- AKANTU_DEBUG_IN();
+// /* -------------------------------------------------------------------------- */
+// Array<Real> & operator*=(Array<Real> & vect, const PETScMatrix & mat) {
+// AKANTU_DEBUG_IN();
- Array<Real> result
- MatMult(this->mat, vect, )
+// Array<Real> result
+// MatMult(this->mat, vect, )
- // AKANTU_DEBUG_ASSERT((vect.getSize()*vect.getNbComponent() == mat.getSize()) &&
- // (vect.getNbComponent() == mat.getNbDegreOfFreedom()),
- // "The size of the matrix and the vector do not match");
+// // AKANTU_DEBUG_ASSERT((vect.getSize()*vect.getNbComponent() == mat.getSize()) &&
+// // (vect.getNbComponent() == mat.getNbDegreOfFreedom()),
+// // "The size of the matrix and the vector do not match");
- const PetscMatrixType & sparse_matrix_type = mat.getPetscMatrixType();
- DOFSynchronizer * dof_synchronizer = mat.getDOFSynchronizerPointer();
+// const PETScMatrixType & sparse_matrix_type = mat.getPETScMatrixType();
+// DOFSynchronizer * dof_synchronizer = mat.getDOFSynchronizerPointer();
- UInt nb_non_zero = mat.getNbNonZero();
- Real * tmp = new Real [vect.getNbComponent() * vect.getSize()];
- std::fill_n(tmp, vect.getNbComponent() * vect.getSize(), 0);
+// UInt nb_non_zero = mat.getNbNonZero();
+// Real * tmp = new Real [vect.getNbComponent() * vect.getSize()];
+// std::fill_n(tmp, vect.getNbComponent() * vect.getSize(), 0);
- Int * i_val = mat.getIRN().storage();
- Int * j_val = mat.getJCN().storage();
- Real * a_val = mat.getA().storage();
+// Int * i_val = mat.getIRN().storage();
+// Int * j_val = mat.getJCN().storage();
+// Real * a_val = mat.getA().storage();
- Real * vect_val = vect.storage();
+// Real * vect_val = vect.storage();
- for (UInt k = 0; k < nb_non_zero; ++k) {
- UInt i = *(i_val++);
- UInt j = *(j_val++);
- Real a = *(a_val++);
+// for (UInt k = 0; k < nb_non_zero; ++k) {
+// UInt i = *(i_val++);
+// UInt j = *(j_val++);
+// Real a = *(a_val++);
- UInt local_i = i - 1;
- UInt local_j = j - 1;
- if(dof_synchronizer) {
- local_i = dof_synchronizer->getDOFLocalID(local_i);
- local_j = dof_synchronizer->getDOFLocalID(local_j);
- }
+// UInt local_i = i - 1;
+// UInt local_j = j - 1;
+// if(dof_synchronizer) {
+// local_i = dof_synchronizer->getDOFLocalID(local_i);
+// local_j = dof_synchronizer->getDOFLocalID(local_j);
+// }
- tmp[local_i] += a * vect_val[local_j];
- if((sparse_matrix_type == _symmetric) && (local_i != local_j))
- tmp[local_j] += a * vect_val[local_i];
- }
+// tmp[local_i] += a * vect_val[local_j];
+// if((sparse_matrix_type == _symmetric) && (local_i != local_j))
+// tmp[local_j] += a * vect_val[local_i];
+// }
- memcpy(vect_val, tmp, vect.getNbComponent() * vect.getSize() * sizeof(Real));
- delete [] tmp;
+// memcpy(vect_val, tmp, vect.getNbComponent() * vect.getSize() * sizeof(Real));
+// delete [] tmp;
- if(dof_synchronizer)
- dof_synchronizer->reduceSynchronize<AddOperation>(vect);
+// if(dof_synchronizer)
+// dof_synchronizer->reduceSynchronize<AddOperation>(vect);
- AKANTU_DEBUG_OUT();
+// AKANTU_DEBUG_OUT();
- return vect;
-}
+// return vect;
+// }
-/* -------------------------------------------------------------------------- */
-void PetscMatrix::copyContent(const PetscMatrix & matrix) {
- AKANTU_DEBUG_IN();
- AKANTU_DEBUG_ASSERT(nb_non_zero == matrix.getNbNonZero(),
- "The two matrices don't have the same profiles");
+// /* -------------------------------------------------------------------------- */
+// void PETScMatrix::copyContent(const PETScMatrix & matrix) {
+// AKANTU_DEBUG_IN();
+// AKANTU_DEBUG_ASSERT(nb_non_zero == matrix.getNbNonZero(),
+// "The two matrices don't have the same profiles");
- MatCopy(this->mat, matrix->mat, SAME_NONZERO_PATTERN);
+// MatCopy(this->mat, matrix->mat, SAME_NONZERO_PATTERN);
- AKANTU_DEBUG_OUT();
-}
+// AKANTU_DEBUG_OUT();
+// }
-///* -------------------------------------------------------------------------- */
-//void PetscMatrix::copyProfile(const PetscMatrix & matrix) {
-// AKANTU_DEBUG_IN();
-// irn = matrix.irn;
-// jcn = matrix.jcn;
-// nb_non_zero = matrix.nb_non_zero;
-// irn_jcn_k = matrix.irn_jcn_k;
-// a.resize(nb_non_zero);
-// AKANTU_DEBUG_OUT();
-//}
+// ///* -------------------------------------------------------------------------- */
+// //void PETScMatrix::copyProfile(const PETScMatrix & matrix) {
+// // AKANTU_DEBUG_IN();
+// // irn = matrix.irn;
+// // jcn = matrix.jcn;
+// // nb_non_zero = matrix.nb_non_zero;
+// // irn_jcn_k = matrix.irn_jcn_k;
+// // a.resize(nb_non_zero);
+// // AKANTU_DEBUG_OUT();
+// //}
/* -------------------------------------------------------------------------- */
-void PetscMatrix::add(const SparseMatrix & matrix, Real alpha) {
+void PETScMatrix::add(const SparseMatrix & matrix, Real alpha) {
+ PetscErrorCode ierr;
AKANTU_DEBUG_ASSERT(nb_non_zero == matrix.getNbNonZero(),
- "The two matrices don't have the same profile");
-
- // get values and local indices of matrix that has to be added
- UInt i_global;
- UInt j_global;
+ "The two matrix don't have the same profiles");
- Array<UInt>::scalar_iterator i_global_it = i_global.begin();
- Array<UInt>::scalar_iterator i_global_end = i_global.end();
- Array<UInt>::scalar_iterator j_global_it = j_global.begin();
+ const Array<Real> matrix_a = matrix.getA();
+ Array<Int> matrix_irn = matrix.getIRN();
+ Array<Int> matrix_jcn = matrix.getJCN();
+ UInt matrix_nb_non_zero = matrix.getNbNonZero();
- Int * i = matrix.irn.storage();
- Int *j = matrix.jcn.storage();
- Real * a_val = matrix.a.storage();
+ /// map global akantu to global petsc
+ ierr = AOApplicationToPetsc(this->petsc_wrapper->ao, matrix_nb_non_zero,matrix_irn.storage()); CHKERRXX(ierr);
+ ierr = AOApplicationToPetsc(this->petsc_wrapper->ao, matrix_nb_non_zero,matrix_jcn.storage()); CHKERRXX(ierr);
- for (; i != i_global_end; ++i_global_it; ++j_global_it; ++i; ++j; ++a_val) {
- // get the akantu global dof index
- i_global = dof_synchronizer->getDOFGlobalID(*i);
- j_global = dof_synchronizer->getDOFGlobalID(*j);
+ Real val_to_add = 0;
- // get the petsc gloobal dof index
- AOApplicationToPetsc(&ao,2,{i_global,j_global})
- MatSetValue(*mat,i,j,*a_val,INSERT_VALUES);
+ for (UInt n = 0; n < matrix_nb_non_zero; ++n) {
+ val_to_add = alpha * matrix_a(n);
+ ierr = MatSetValue(this->petsc_wrapper->mat, matrix_irn(n), matrix_jcn(n), val_to_add,INSERT_VALUES); CHKERRXX(ierr);
}
-
-}
+}
/* -------------------------------------------------------------------------- */
-void PetscMatrix::lump(Array<Real> & lumped) {
- AKANTU_DEBUG_IN();
-
- AKANTU_DEBUG_ASSERT((lumped.getNbComponent() == nb_degree_of_freedom),
- "The size of the matrix and the vector do not match");
-
- UInt vect_size = size / nb_degree_of_freedom;
- if(dof_synchronizer) vect_size = dof_synchronizer->getNbDOFs() / nb_degree_of_freedom;
-
- lumped.resize(vect_size);
- lumped.clear();
-
- Int * i_val = irn.storage();
- Int * j_val = jcn.storage();
- Real * a_val = a.storage();
-
- Real * vect_val = lumped.storage();
-
- for (UInt k = 0; k < nb_non_zero; ++k) {
- UInt i = *(i_val++);
- UInt j = *(j_val++);
- Real a = *(a_val++);
-
- UInt local_i = i - 1;
- UInt local_j = j - 1;
- if(dof_synchronizer) {
- local_i = dof_synchronizer->getDOFLocalID(local_i);
- local_j = dof_synchronizer->getDOFLocalID(local_j);
- }
-
- vect_val[local_i] += a;
- if(sparse_matrix_type == _symmetric && (i != j))
- vect_val[local_j] += a;
- }
-
- if(dof_synchronizer)
- dof_synchronizer->reduceSynchronize<AddOperation>(lumped);
-
- AKANTU_DEBUG_OUT();
+void PETScMatrix::performAssembly() {
+ PetscErrorCode ierr;
+ ierr = MatAssemblyBegin(this->petsc_wrapper->mat, MAT_FINAL_ASSEMBLY); CHKERRXX(ierr);
+ ierr = MatAssemblyEnd(this->petsc_wrapper->mat, MAT_FINAL_ASSEMBLY); CHKERRXX(ierr);
}
+
__END_AKANTU__
diff --git a/src/solver/petsc_matrix.hh b/src/solver/petsc_matrix.hh
index 75d7bc308..40a04103d 100644
--- a/src/solver/petsc_matrix.hh
+++ b/src/solver/petsc_matrix.hh
@@ -1,126 +1,158 @@
/**
* @file petsc_matrix.hh
* @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @date Mon Jul 21 14:49:49 2014
*
* @brief Interface for PETSc matrices
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_PETSC_MATRIX_HH__
#define __AKANTU_PETSC_MATRIX_HH__
/* -------------------------------------------------------------------------- */
#include "sparse_matrix.hh"
-struct Mat;
-struct AO;
-
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
+class PETScWrapper;
+
+class PETScMatrix : public SparseMatrix {
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+public:
+ PETScMatrix(UInt size,
+ const SparseMatrixType & sparse_matrix_type,
+ const ID & id = "petsc_matrix",
+ const MemoryID & memory_id = 0);
+
+ PETScMatrix(const SparseMatrix & matrix,
+ const ID & id = "petsc_matrix",
+ const MemoryID & memory_id = 0);
+ virtual ~PETScMatrix();
+
+
+private:
+ /// init internal PETSc matrix
+ void init();
-class PetscMatrix : SparseMatrix {
/* ------------------------------------------------------------------------ */
- /* Typedefs */
- /* ------------------------------------------------------------------------ */
+ /* Methods */
+ /* ------------------------------------------------------------------------ */
+public:
+ /// resize PETSc matrix
+ void resize(const DOFSynchronizer & dof_synchronizer);
- class PetscMatrixValue {
+ // /// remove the existing profile
+ // inline void clearProfile();
- public:
+ // /// add a non-zero element
+ // virtual UInt addToProfile(UInt i, UInt j);
- PetscMatrixValue(){}
+ // /// set the matrix to 0
+ // inline void clear();
- inline PetscMatrixValue operator = (Real v){
- MatSetValue(*mat,i,j,v,INSERT_VALUES);
- return *this;
- };
+ // /// assemble a local matrix in the sparse one
+ // inline void addToMatrix(UInt i, UInt j, Real value);
+ /// fill the profil of the matrix
+ virtual void buildProfile(const Mesh & mesh, const DOFSynchronizer & dof_synchronizer, UInt nb_degree_of_freedom);
- inline PetscMatrixValue operator += (Real v){
- // std::cout << "set " << std::abs(Real(i)-Real(j)) << std::endl;
- MatSetValue(*mat,i,j,v,ADD_VALUES);
- return *this;
- };
+ // /// modify the matrix to "remove" the blocked dof
+ // virtual void applyBoundary(const Array<bool> & boundary, Real block_val = 1.);
- inline PetscMatrixValue operator -= (Real v){
- MatSetValue(*mat,i,j,-1.*v,ADD_VALUES);
- return *this;
- };
+ /// modify the matrix to "remove" the blocked dof
+ virtual void applyBoundaryNormal(Array<bool> & boundary_normal, Array<Real> & EulerAngles, Array<Real> & rhs, const Array<Real> & matrix, Array<Real> & rhs_rotated) {
+ AKANTU_DEBUG_TO_IMPLEMENT();
+ }
- friend PetscMatrix;
+ /// modify the matrix to "remove" the blocked dof
+ virtual void removeBoundary(const Array<bool> & boundary) {
+ AKANTU_DEBUG_TO_IMPLEMENT();
+ }
+
+ /// perform assembly so that matrix is ready for use
+ void performAssembly();
- private:
+ // /// restore the profile that was before removing the boundaries
+ // virtual void restoreProfile();
- UInt i,j;
- Mat * mat;
- };
+ // /// save the profil in a file using the MatrixMarket file format
+ // virtual void saveProfile(const std::string & filename) const;
- /* ------------------------------------------------------------------------ */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------ */
-public:
- PetscMatrix(UInt size,
- const SparseMatrixType & sparse_matrix_type,
- UInt nb_degree_of_freedom,
- const ID & id = "petsc_matrix",
- const MemoryID & memory_id = 0);
+ /// save the matrix in a file using the MatrixMarket file format
+ virtual void saveMatrix(const std::string & filename) const;
- PetscMatrix(const PetscMatrix & matrix,
- const ID & id = "petsc_matrix",
- const MemoryID & memory_id = 0);
+ // /// copy assuming the profile are the same
+ // virtual void copyContent(const SparseMatrix & matrix);
- virtual ~PetscMatrix();
+ // /// copy profile
+ // // void copyProfile(const SparseMatrix & matrix);
- typedef std::pair<UInt, UInt> KeyCOO;
- typedef unordered_map<KeyCOO, UInt>::type coordinate_list_map;
+ /// add matrix assuming the profile are the same
+ virtual void add(const SparseMatrix & matrix, Real alpha);
+ // /// diagonal lumping
+ // virtual void lump(Array<Real> & lumped);
- /* ------------------------------------------------------------------------ */
- /* Methods */
- /* ------------------------------------------------------------------------ */
-public:
+ // /// function to print the contain of the class
+ // //virtual void printself(std::ostream & stream, int indent = 0) const;
+
+private:
+ /// create a mapping from the global akantu numbering to the PETSc numbering
+ void createGlobalAkantuToPETScMap(Int* local_master_eq_nbs_ptr);
+
+ /// create a mapping from the local akantu numbering to the PETSc numbering
+ void createLocalAkantuToPETScMap(const DOFSynchronizer & dof_synchronizer);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
-
+
private:
+ /// store the PETSc structures
+ PETScWrapper * petsc_wrapper;
- Mat * mat;
- AO * ao;
+ /// size of the diagonal part of the matrix partition
+ Int local_size;
+ /// number of nonzeros in every row of the diagonal part
+ Array<Int> d_nnz;
+ /// number of nonzeros in every row of the off-diagonal part
+ Array<Int> o_nnz;
};
__END_AKANTU__
#endif /* __AKANTU_PETSCI_MATRIX_HH__ */
diff --git a/src/solver/petsc_matrix_inline_impl.cc b/src/solver/petsc_matrix_inline_impl.cc
index a4c0792ae..b3fe55d54 100644
--- a/src/solver/petsc_matrix_inline_impl.cc
+++ b/src/solver/petsc_matrix_inline_impl.cc
@@ -1,73 +1,73 @@
/**
* @file petsc_matrix_inline_impl.cc
* @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @date Mon Jul 21 21:42:12 2014
*
* @brief Implementation of PetscMatrix inline functions
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
-/* -------------------------------------------------------------------------- */
-inline UInt PetscMatrix::addToProfile(UInt i, UInt j) {
+// /* -------------------------------------------------------------------------- */
+// inline UInt PetscMatrix::addToProfile(UInt i, UInt j) {
- // initialize new value
- MatSetValuesLocal(this->mat,1,i,1,j,0,INSERT_VALUES);
+// // initialize new value
+// MatSetValuesLocal(this->mat,1,i,1,j,0,INSERT_VALUES);
- nb_non_zero++;
+// nb_non_zero++;
- return nb_non_zero - 1;
-}
+// return nb_non_zero - 1;
+// }
-/* -------------------------------------------------------------------------- */
-inline void PetscMatrix::clearProfile() {
- //retain matrix but reinitialize its content
- MatZeroEntries(this->mat);
- nb_non_zero = 0;
-}
+// /* -------------------------------------------------------------------------- */
+// inline void PetscMatrix::clearProfile() {
+// //retain matrix but reinitialize its content
+// MatZeroEntries(this->mat);
+// nb_non_zero = 0;
+// }
-/* -------------------------------------------------------------------------- */
-inline void PetscMatrix::clear() {
- memset(a.storage(), 0, nb_non_zero*sizeof(Real));
-}
+// /* -------------------------------------------------------------------------- */
+// inline void PetscMatrix::clear() {
+// memset(a.storage(), 0, nb_non_zero*sizeof(Real));
+// }
-/* -------------------------------------------------------------------------- */
-inline void PetscMatrix::addToMatrix(UInt i, UInt j, Real value) {
- MatSetValue(this->mat,1,i,1,j,value,ADD_VALUES);
-}
+// /* -------------------------------------------------------------------------- */
+// inline void PetscMatrix::addToMatrix(UInt i, UInt j, Real value) {
+// MatSetValue(this->mat,1,i,1,j,value,ADD_VALUES);
+// }
-/* -------------------------------------------------------------------------- */
-inline Real PetscMatrix::operator()(UInt i, UInt j) const {
- KeyCOO jcn_irn = key(i, j);
- coordinate_list_map::const_iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn);
- if(irn_jcn_k_it == irn_jcn_k.end()) return 0;
- return this->mat(i,j).
-}
+// /* -------------------------------------------------------------------------- */
+// inline Real PetscMatrix::operator()(UInt i, UInt j) const {
+// KeyCOO jcn_irn = key(i, j);
+// coordinate_list_map::const_iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn);
+// if(irn_jcn_k_it == irn_jcn_k.end()) return 0;
+// return this->mat(i,j).
+// }
-/* -------------------------------------------------------------------------- */
-inline Real & PetscMatrix::operator()(UInt i, UInt j) {
- KeyCOO jcn_irn = key(i, j);
- coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn);
- AKANTU_DEBUG_ASSERT(irn_jcn_k_it != irn_jcn_k.end(),
- "Couple (i,j) = (" << i << "," << j << ") does not exist in the profile");
+// /* -------------------------------------------------------------------------- */
+// inline Real & PetscMatrix::operator()(UInt i, UInt j) {
+// KeyCOO jcn_irn = key(i, j);
+// coordinate_list_map::iterator irn_jcn_k_it = irn_jcn_k.find(jcn_irn);
+// AKANTU_DEBUG_ASSERT(irn_jcn_k_it != irn_jcn_k.end(),
+// "Couple (i,j) = (" << i << "," << j << ") does not exist in the profile");
- return a.storage()[irn_jcn_k_it->second];
-}
+// return a.storage()[irn_jcn_k_it->second];
+// }
diff --git a/src/solver/solver_petsc.cc b/src/solver/solver_petsc.cc
index a2bd2348e..ce058bc87 100644
--- a/src/solver/solver_petsc.cc
+++ b/src/solver/solver_petsc.cc
@@ -1,954 +1,954 @@
-/**
- * @file solver_petsc.cc
- *
- * @author Nicolas Richart <nicolas.richart@epfl.ch>
- # @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
- *
- * @date Mon Dec 13 10:48:06 2010
- *
- * @brief Solver class implementation for the petsc solver
- *
- * @section LICENSE
- *
- * Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
- *
- */
+// /**
+// * @file solver_petsc.cc
+// *
+// * @author Nicolas Richart <nicolas.richart@epfl.ch>
+// # @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
+// *
+// * @date Mon Dec 13 10:48:06 2010
+// *
+// * @brief Solver class implementation for the petsc solver
+// *
+// * @section LICENSE
+// *
+// * Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
+// *
+// */
-/* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_MPI
-#include "static_communicator_mpi.hh"
-#endif
+// /* -------------------------------------------------------------------------- */
+// #ifdef AKANTU_USE_MPI
+// #include "static_communicator_mpi.hh"
+// #endif
-#include "solver_petsc.hh"
-#include "dof_synchronizer.hh"
+// #include "solver_petsc.hh"
+// #include "dof_synchronizer.hh"
-__BEGIN_AKANTU__
+// __BEGIN_AKANTU__
-void finalize_petsc() {
+// void finalize_petsc() {
- static bool finalized = false;
- if (!finalized) {
+// static bool finalized = false;
+// if (!finalized) {
- cout<<"*** INFO *** PETSc is finalizing..."<<endl;
- // finalize PETSc
- PetscErrorCode ierr = PetscFinalize();CHKERRCONTINUE(ierr);
- finalized = true;
- }
-}
+// cout<<"*** INFO *** PETSc is finalizing..."<<endl;
+// // finalize PETSc
+// PetscErrorCode ierr = PetscFinalize();CHKERRCONTINUE(ierr);
+// finalized = true;
+// }
+// }
-SolverPETSc::sparse_vector_type
-SolverPETSc::operator()(const SolverPETSc::sparse_matrix_type& AA,
- const SolverPETSc::sparse_vector_type& bb) {
-
-
-#ifdef CPPUTILS_VERBOSE
- // parallel output stream
- Output_stream out;
- // timer
- cpputils::ctimer timer;
- out<<"Inside PETSc solver: "<<timer<<endl;
-#endif
-
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Inside operator()(const sparse_matrix_type&, const sparse_vector_type&)... "<<timer<<endl;
-#endif
-
- assert(AA.rows() == bb.size());
-
- // KSP ksp; //!< linear solver context
-
- Vec b; /* RHS */
- PC pc; /* preconditioner context */
- PetscErrorCode ierr;
- PetscInt nlocal;
- PetscInt n = bb.size();
- VecScatter ctx;
-
- /*
- Create vectors. Note that we form 1 vector from scratch and
- then duplicate as needed. For this simple case let PETSc decide how
- many elements of the vector are stored on each processor. The second
- argument to VecSetSizes() below causes PETSc to decide.
- */
- ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRCONTINUE(ierr);
- ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRCONTINUE(ierr);
- ierr = VecSetFromOptions(b);CHKERRCONTINUE(ierr);
- if (!allocated_) {
- ierr = VecDuplicate(b,&x_);CHKERRCONTINUE(ierr);
- } else
- VecZeroEntries(x_);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Vectors created... "<<timer<<endl;
-#endif
-
-
- /* Set hight-hand-side vector */
- for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
- int row = it->first;
- ierr = VecSetValues(b, 1, &row, &it->second, ADD_VALUES);
- }
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Right hand side set... "<<timer<<endl;
-#endif
-
-
- /*
- Assemble vector, using the 2-step process:
- VecAssemblyBegin(), VecAssemblyEnd()
- Computations can be done while messages are in transition
- by placing code between these two statements.
- */
- ierr = VecAssemblyBegin(b);CHKERRCONTINUE(ierr);
- ierr = VecAssemblyEnd(b);CHKERRCONTINUE(ierr);
-
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Right-hand-side vector assembled... "<<timer<<endl;
-
- ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRCONTINUE(ierr);
-
- Vec b_all;
- ierr = VecScatterCreateToAll(b, &ctx, &b_all);CHKERRCONTINUE(ierr);
- ierr = VecScatterBegin(ctx,b,b_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
- ierr = VecScatterEnd(ctx,b,b_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
-
- value_type nrm;
- VecNorm(b_all,NORM_2,&nrm);
- out<<" Norm of right hand side... "<<nrm<<endl;
-#endif
-
- /* Identify the starting and ending mesh points on each
- processor for the interior part of the mesh. We let PETSc decide
- above. */
-
- // PetscInt rstart,rend;
- // ierr = VecGetOwnershipRange(x_,&rstart,&rend);CHKERRCONTINUE(ierr);
- ierr = VecGetLocalSize(x_,&nlocal);CHKERRCONTINUE(ierr);
-
-
- /*
- Create matrix. When using MatCreate(), the matrix format can
- be specified at runtime.
+// SolverPETSc::sparse_vector_type
+// SolverPETSc::operator()(const SolverPETSc::sparse_matrix_type& AA,
+// const SolverPETSc::sparse_vector_type& bb) {
+
+
+// #ifdef CPPUTILS_VERBOSE
+// // parallel output stream
+// Output_stream out;
+// // timer
+// cpputils::ctimer timer;
+// out<<"Inside PETSc solver: "<<timer<<endl;
+// #endif
+
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Inside operator()(const sparse_matrix_type&, const sparse_vector_type&)... "<<timer<<endl;
+// #endif
+
+// assert(AA.rows() == bb.size());
+
+// // KSP ksp; //!< linear solver context
+
+// Vec b; /* RHS */
+// PC pc; /* preconditioner context */
+// PetscErrorCode ierr;
+// PetscInt nlocal;
+// PetscInt n = bb.size();
+// VecScatter ctx;
+
+// /*
+// Create vectors. Note that we form 1 vector from scratch and
+// then duplicate as needed. For this simple case let PETSc decide how
+// many elements of the vector are stored on each processor. The second
+// argument to VecSetSizes() below causes PETSc to decide.
+// */
+// ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRCONTINUE(ierr);
+// ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRCONTINUE(ierr);
+// ierr = VecSetFromOptions(b);CHKERRCONTINUE(ierr);
+// if (!allocated_) {
+// ierr = VecDuplicate(b,&x_);CHKERRCONTINUE(ierr);
+// } else
+// VecZeroEntries(x_);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Vectors created... "<<timer<<endl;
+// #endif
+
+
+// /* Set hight-hand-side vector */
+// for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
+// int row = it->first;
+// ierr = VecSetValues(b, 1, &row, &it->second, ADD_VALUES);
+// }
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Right hand side set... "<<timer<<endl;
+// #endif
+
+
+// /*
+// Assemble vector, using the 2-step process:
+// VecAssemblyBegin(), VecAssemblyEnd()
+// Computations can be done while messages are in transition
+// by placing code between these two statements.
+// */
+// ierr = VecAssemblyBegin(b);CHKERRCONTINUE(ierr);
+// ierr = VecAssemblyEnd(b);CHKERRCONTINUE(ierr);
+
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Right-hand-side vector assembled... "<<timer<<endl;
+
+// ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRCONTINUE(ierr);
+
+// Vec b_all;
+// ierr = VecScatterCreateToAll(b, &ctx, &b_all);CHKERRCONTINUE(ierr);
+// ierr = VecScatterBegin(ctx,b,b_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
+// ierr = VecScatterEnd(ctx,b,b_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
+
+// value_type nrm;
+// VecNorm(b_all,NORM_2,&nrm);
+// out<<" Norm of right hand side... "<<nrm<<endl;
+// #endif
+
+// /* Identify the starting and ending mesh points on each
+// processor for the interior part of the mesh. We let PETSc decide
+// above. */
+
+// // PetscInt rstart,rend;
+// // ierr = VecGetOwnershipRange(x_,&rstart,&rend);CHKERRCONTINUE(ierr);
+// ierr = VecGetLocalSize(x_,&nlocal);CHKERRCONTINUE(ierr);
+
+
+// /*
+// Create matrix. When using MatCreate(), the matrix format can
+// be specified at runtime.
- Performance tuning note: For problems of substantial size,
- preallocation of matrix memory is crucial for attaining good
- performance. See the matrix chapter of the users manual for details.
+// Performance tuning note: For problems of substantial size,
+// preallocation of matrix memory is crucial for attaining good
+// performance. See the matrix chapter of the users manual for details.
- We pass in nlocal as the "local" size of the matrix to force it
- to have the same parallel layout as the vector created above.
- */
- if (!allocated_) {
+// We pass in nlocal as the "local" size of the matrix to force it
+// to have the same parallel layout as the vector created above.
+// */
+// if (!allocated_) {
- ierr = MatCreate(PETSC_COMM_WORLD,&A_);CHKERRCONTINUE(ierr);
- ierr = MatSetSizes(A_,nlocal,nlocal,n,n);CHKERRCONTINUE(ierr);
- ierr = MatSetFromOptions(A_);CHKERRCONTINUE(ierr);
- ierr = MatSetUp(A_);CHKERRCONTINUE(ierr);
- } else {
- // zero-out matrix
- MatZeroEntries(A_);
- }
-
-
- /*
- Assemble matrix.
+// ierr = MatCreate(PETSC_COMM_WORLD,&A_);CHKERRCONTINUE(ierr);
+// ierr = MatSetSizes(A_,nlocal,nlocal,n,n);CHKERRCONTINUE(ierr);
+// ierr = MatSetFromOptions(A_);CHKERRCONTINUE(ierr);
+// ierr = MatSetUp(A_);CHKERRCONTINUE(ierr);
+// } else {
+// // zero-out matrix
+// MatZeroEntries(A_);
+// }
+
+
+// /*
+// Assemble matrix.
- The linear system is distributed across the processors by
- chunks of contiguous rows, which correspond to contiguous
- sections of the mesh on which the problem is discretized.
- For matrix assembly, each processor contributes entries for
- the part that it owns locally.
- */
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Zeroed-out sparse matrix entries... "<<timer<<endl;
-#endif
-
- for (sparse_matrix_type::const_hash_iterator it = AA.map_.begin(); it != AA.map_.end(); ++it) {
+// The linear system is distributed across the processors by
+// chunks of contiguous rows, which correspond to contiguous
+// sections of the mesh on which the problem is discretized.
+// For matrix assembly, each processor contributes entries for
+// the part that it owns locally.
+// */
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Zeroed-out sparse matrix entries... "<<timer<<endl;
+// #endif
+
+// for (sparse_matrix_type::const_hash_iterator it = AA.map_.begin(); it != AA.map_.end(); ++it) {
- // get subscripts
- std::pair<size_t,size_t> subs = AA.unhash(it->first);
- PetscInt row = subs.first;
- PetscInt col = subs.second;
- ierr = MatSetValues(A_, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERRCONTINUE(ierr);
- }
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Filled sparse matrix... "<<timer<<endl;
-#endif
-
-
-
-
- /* Assemble the matrix */
- ierr = MatAssemblyBegin(A_,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
- ierr = MatAssemblyEnd(A_,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
-
-
- if (!allocated_) {
- // set after the first MatAssemblyEnd
- // ierr = MatSetOption(A_, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);CHKERRCONTINUE(ierr);
- ierr = MatSetOption(A_, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRCONTINUE(ierr);
- }
-
-
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Sparse matrix assembled... "<<timer<<endl;
- // view matrix
- MatView(A_, PETSC_VIEWER_STDOUT_WORLD);
-
- MatNorm(A_,NORM_FROBENIUS,&nrm);
- out<<" Norm of stiffness matrix... "<<nrm<<endl;
-#endif
-
- /*
- Set operators. Here the matrix that defines the linear system
- also serves as the preconditioning matrix.
- */
- // ierr = KSPSetOperators(ksp,A_,A_,DIFFERENT_NONZERO_PATTERN);CHKERRCONTINUE(ierr);
- ierr = KSPSetOperators(ksp_,A_,A_,SAME_NONZERO_PATTERN);CHKERRCONTINUE(ierr);
-
-
- //
- // /*
- // Set runtime options, e.g.,
- // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
- // These options will override those specified above as long as
- // KSPSetFromOptions() is called _after_ any other customization
- // routines.
- // */
- // ierr = KSPSetFromOptions(ksp);CHKERRCONTINUE(ierr);
-
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Solving system... "<<timer<<endl;
-#endif
-
-
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
- /*
- Solve linear system
- */
- ierr = KSPSolve(ksp_,b,x_);CHKERRCONTINUE(ierr);
-
-
-#ifdef CPPUTILS_VERBOSE
-
- /*
- View solver info; we could instead use the option -ksp_view to
- print this info to the screen at the conclusion of KSPSolve().
- */
- ierr = KSPView(ksp_,PETSC_VIEWER_STDOUT_WORLD);CHKERRCONTINUE(ierr);
-
- int iter;
- KSPGetIterationNumber(ksp_, &iter);
- out<<" System solved in "<<iter<<" iterations... "<<timer<<endl;
- KSPConvergedReason reason;
- ierr = KSPGetConvergedReason(ksp_,&reason);CHKERRCONTINUE(ierr);
- if (reason < 0)
- out<<"*** WARNING *** PETSc solver diverged with flag ";
- else
- out<<"*** INFO *** PETSc solver converged with flag ";
-
- if (reason == KSP_CONVERGED_RTOL)
- out<<"KSP_CONVERGED_RTOL"<<endl;
- else if (reason == KSP_CONVERGED_ATOL)
- out<<"KSP_CONVERGED_ATOL"<<endl;
- else if (reason == KSP_CONVERGED_ITS)
- out<<"KSP_CONVERGED_ITS"<<endl;
- else if (reason == KSP_CONVERGED_CG_NEG_CURVE)
- out<<"KSP_CONVERGED_CG_NEG_CURVE"<<endl;
- else if (reason == KSP_CONVERGED_CG_CONSTRAINED)
- out<<"KSP_CONVERGED_CG_CONSTRAINED"<<endl;
- else if (reason == KSP_CONVERGED_STEP_LENGTH)
- out<<"KSP_CONVERGED_STEP_LENGTH"<<endl;
- else if (reason == KSP_CONVERGED_HAPPY_BREAKDOWN)
- out<<"KSP_CONVERGED_HAPPY_BREAKDOWN"<<endl;
- else if (reason == KSP_DIVERGED_NULL)
- out<<"KSP_DIVERGED_NULL"<<endl;
- else if (reason == KSP_DIVERGED_ITS)
- out<<"KSP_DIVERGED_ITS"<<endl;
- else if (reason == KSP_DIVERGED_DTOL)
- out<<"KSP_DIVERGED_DTOL"<<endl;
- else if (reason == KSP_DIVERGED_BREAKDOWN)
- out<<"KSP_DIVERGED_BREAKDOWN"<<endl;
- else if (reason == KSP_DIVERGED_BREAKDOWN_BICG)
- out<<"KSP_DIVERGED_BREAKDOWN_BICG"<<endl;
- else if (reason == KSP_DIVERGED_NONSYMMETRIC)
- out<<"KSP_DIVERGED_NONSYMMETRIC"<<endl;
- else if (reason == KSP_DIVERGED_INDEFINITE_PC)
- out<<"KSP_DIVERGED_INDEFINITE_PC"<<endl;
- else if (reason == KSP_DIVERGED_NAN)
- out<<"KSP_DIVERGED_NAN"<<endl;
- else if (reason == KSP_DIVERGED_INDEFINITE_MAT)
- out<<"KSP_DIVERGED_INDEFINITE_MAT"<<endl;
- else if (reason == KSP_CONVERGED_ITERATING)
- out<<"KSP_CONVERGED_ITERATING"<<endl;
-
- PetscReal rnorm;
- ierr = KSPGetResidualNorm(ksp_,&rnorm);CHKERRCONTINUE(ierr);
-
- out<<"PETSc last residual norm computed: "<<rnorm<<endl;
-
- ierr = VecView(x_,PETSC_VIEWER_STDOUT_WORLD);CHKERRCONTINUE(ierr);
-
- VecNorm(x_,NORM_2,&nrm);
- out<<" Norm of solution vector... "<<nrm<<endl;
-
-#endif
-
-
-
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Check solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
- Vec x_all;
-
- ierr = VecScatterCreateToAll(x_, &ctx, &x_all);CHKERRCONTINUE(ierr);
- ierr = VecScatterBegin(ctx,x_,x_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
- ierr = VecScatterEnd(ctx,x_,x_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
-
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Solution vector scattered... "<<timer<<endl;
- VecNorm(x_all,NORM_2,&nrm);
- out<<" Norm of scattered vector... "<<nrm<<endl;
- // ierr = VecView(x_all,PETSC_VIEWER_STDOUT_WORLD);CHKERRCONTINUE(ierr);
-#endif
-
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Get values from solution and store them in the object that will be
- returned
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
-
- sparse_vector_type xx(bb.size());
-
- /* Set solution vector */
- double zero = 0.;
- double val;
- for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
- int row = it->first;
- ierr = VecGetValues(x_all, 1, &row, &val);
- if (val != zero)
- xx[row] = val;
- }
-
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Solution vector copied... "<<timer<<endl;
- // out<<" Norm of copied solution vector... "<<norm(xx, Insert_t)<<endl;
-#endif
-
-
- /*
- Free work space. All PETSc objects should be destroyed when they
- are no longer needed.
- */
- ierr = VecDestroy(&b);CHKERRCONTINUE(ierr);
- // ierr = KSPDestroy(&ksp);CHKERRCONTINUE(ierr);
-
- // set allocated flag
- if (!allocated_) {
- allocated_ = true;
- }
-
-
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Temporary data structures destroyed... "<<timer<<endl;
-#endif
-
- return xx;
-}
+// // get subscripts
+// std::pair<size_t,size_t> subs = AA.unhash(it->first);
+// PetscInt row = subs.first;
+// PetscInt col = subs.second;
+// ierr = MatSetValues(A_, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERRCONTINUE(ierr);
+// }
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Filled sparse matrix... "<<timer<<endl;
+// #endif
+
+
+
+
+// /* Assemble the matrix */
+// ierr = MatAssemblyBegin(A_,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
+// ierr = MatAssemblyEnd(A_,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
+
+
+// if (!allocated_) {
+// // set after the first MatAssemblyEnd
+// // ierr = MatSetOption(A_, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);CHKERRCONTINUE(ierr);
+// ierr = MatSetOption(A_, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);CHKERRCONTINUE(ierr);
+// }
+
+
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Sparse matrix assembled... "<<timer<<endl;
+// // view matrix
+// MatView(A_, PETSC_VIEWER_STDOUT_WORLD);
+
+// MatNorm(A_,NORM_FROBENIUS,&nrm);
+// out<<" Norm of stiffness matrix... "<<nrm<<endl;
+// #endif
+
+// /*
+// Set operators. Here the matrix that defines the linear system
+// also serves as the preconditioning matrix.
+// */
+// // ierr = KSPSetOperators(ksp,A_,A_,DIFFERENT_NONZERO_PATTERN);CHKERRCONTINUE(ierr);
+// ierr = KSPSetOperators(ksp_,A_,A_,SAME_NONZERO_PATTERN);CHKERRCONTINUE(ierr);
+
+
+// //
+// // /*
+// // Set runtime options, e.g.,
+// // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
+// // These options will override those specified above as long as
+// // KSPSetFromOptions() is called _after_ any other customization
+// // routines.
+// // */
+// // ierr = KSPSetFromOptions(ksp);CHKERRCONTINUE(ierr);
+
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Solving system... "<<timer<<endl;
+// #endif
+
+
+// /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Solve the linear system
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+// /*
+// Solve linear system
+// */
+// ierr = KSPSolve(ksp_,b,x_);CHKERRCONTINUE(ierr);
+
+
+// #ifdef CPPUTILS_VERBOSE
+
+// /*
+// View solver info; we could instead use the option -ksp_view to
+// print this info to the screen at the conclusion of KSPSolve().
+// */
+// ierr = KSPView(ksp_,PETSC_VIEWER_STDOUT_WORLD);CHKERRCONTINUE(ierr);
+
+// int iter;
+// KSPGetIterationNumber(ksp_, &iter);
+// out<<" System solved in "<<iter<<" iterations... "<<timer<<endl;
+// KSPConvergedReason reason;
+// ierr = KSPGetConvergedReason(ksp_,&reason);CHKERRCONTINUE(ierr);
+// if (reason < 0)
+// out<<"*** WARNING *** PETSc solver diverged with flag ";
+// else
+// out<<"*** INFO *** PETSc solver converged with flag ";
+
+// if (reason == KSP_CONVERGED_RTOL)
+// out<<"KSP_CONVERGED_RTOL"<<endl;
+// else if (reason == KSP_CONVERGED_ATOL)
+// out<<"KSP_CONVERGED_ATOL"<<endl;
+// else if (reason == KSP_CONVERGED_ITS)
+// out<<"KSP_CONVERGED_ITS"<<endl;
+// else if (reason == KSP_CONVERGED_CG_NEG_CURVE)
+// out<<"KSP_CONVERGED_CG_NEG_CURVE"<<endl;
+// else if (reason == KSP_CONVERGED_CG_CONSTRAINED)
+// out<<"KSP_CONVERGED_CG_CONSTRAINED"<<endl;
+// else if (reason == KSP_CONVERGED_STEP_LENGTH)
+// out<<"KSP_CONVERGED_STEP_LENGTH"<<endl;
+// else if (reason == KSP_CONVERGED_HAPPY_BREAKDOWN)
+// out<<"KSP_CONVERGED_HAPPY_BREAKDOWN"<<endl;
+// else if (reason == KSP_DIVERGED_NULL)
+// out<<"KSP_DIVERGED_NULL"<<endl;
+// else if (reason == KSP_DIVERGED_ITS)
+// out<<"KSP_DIVERGED_ITS"<<endl;
+// else if (reason == KSP_DIVERGED_DTOL)
+// out<<"KSP_DIVERGED_DTOL"<<endl;
+// else if (reason == KSP_DIVERGED_BREAKDOWN)
+// out<<"KSP_DIVERGED_BREAKDOWN"<<endl;
+// else if (reason == KSP_DIVERGED_BREAKDOWN_BICG)
+// out<<"KSP_DIVERGED_BREAKDOWN_BICG"<<endl;
+// else if (reason == KSP_DIVERGED_NONSYMMETRIC)
+// out<<"KSP_DIVERGED_NONSYMMETRIC"<<endl;
+// else if (reason == KSP_DIVERGED_INDEFINITE_PC)
+// out<<"KSP_DIVERGED_INDEFINITE_PC"<<endl;
+// else if (reason == KSP_DIVERGED_NAN)
+// out<<"KSP_DIVERGED_NAN"<<endl;
+// else if (reason == KSP_DIVERGED_INDEFINITE_MAT)
+// out<<"KSP_DIVERGED_INDEFINITE_MAT"<<endl;
+// else if (reason == KSP_CONVERGED_ITERATING)
+// out<<"KSP_CONVERGED_ITERATING"<<endl;
+
+// PetscReal rnorm;
+// ierr = KSPGetResidualNorm(ksp_,&rnorm);CHKERRCONTINUE(ierr);
+
+// out<<"PETSc last residual norm computed: "<<rnorm<<endl;
+
+// ierr = VecView(x_,PETSC_VIEWER_STDOUT_WORLD);CHKERRCONTINUE(ierr);
+
+// VecNorm(x_,NORM_2,&nrm);
+// out<<" Norm of solution vector... "<<nrm<<endl;
+
+// #endif
+
+
+
+// /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Check solution and clean up
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+// Vec x_all;
+
+// ierr = VecScatterCreateToAll(x_, &ctx, &x_all);CHKERRCONTINUE(ierr);
+// ierr = VecScatterBegin(ctx,x_,x_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
+// ierr = VecScatterEnd(ctx,x_,x_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
+
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Solution vector scattered... "<<timer<<endl;
+// VecNorm(x_all,NORM_2,&nrm);
+// out<<" Norm of scattered vector... "<<nrm<<endl;
+// // ierr = VecView(x_all,PETSC_VIEWER_STDOUT_WORLD);CHKERRCONTINUE(ierr);
+// #endif
+
+// /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Get values from solution and store them in the object that will be
+// returned
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+// sparse_vector_type xx(bb.size());
+
+// /* Set solution vector */
+// double zero = 0.;
+// double val;
+// for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
+// int row = it->first;
+// ierr = VecGetValues(x_all, 1, &row, &val);
+// if (val != zero)
+// xx[row] = val;
+// }
+
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Solution vector copied... "<<timer<<endl;
+// // out<<" Norm of copied solution vector... "<<norm(xx, Insert_t)<<endl;
+// #endif
+
+
+// /*
+// Free work space. All PETSc objects should be destroyed when they
+// are no longer needed.
+// */
+// ierr = VecDestroy(&b);CHKERRCONTINUE(ierr);
+// // ierr = KSPDestroy(&ksp);CHKERRCONTINUE(ierr);
+
+// // set allocated flag
+// if (!allocated_) {
+// allocated_ = true;
+// }
+
+
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Temporary data structures destroyed... "<<timer<<endl;
+// #endif
+
+// return xx;
+// }
-SolverPETSc::sparse_vector_type SolverPETSc::operator()(const SolverPETSc::sparse_matrix_type& KKpf, const SolverPETSc::sparse_matrix_type& KKpp, const SolverPETSc::sparse_vector_type& UUp) {
+// SolverPETSc::sparse_vector_type SolverPETSc::operator()(const SolverPETSc::sparse_matrix_type& KKpf, const SolverPETSc::sparse_matrix_type& KKpp, const SolverPETSc::sparse_vector_type& UUp) {
-#ifdef CPPUTILS_VERBOSE
- // parallel output stream
- Output_stream out;
- // timer
- cpputils::ctimer timer;
- out<<"Inside SolverPETSc::operator()(const sparse_matrix&, const sparse_matrix&, const sparse_vector&). "<<timer<<endl;
-#endif
+// #ifdef CPPUTILS_VERBOSE
+// // parallel output stream
+// Output_stream out;
+// // timer
+// cpputils::ctimer timer;
+// out<<"Inside SolverPETSc::operator()(const sparse_matrix&, const sparse_matrix&, const sparse_vector&). "<<timer<<endl;
+// #endif
- Mat Kpf, Kpp;
- Vec Up, Pf, Pp;
+// Mat Kpf, Kpp;
+// Vec Up, Pf, Pp;
- PetscErrorCode ierr = MatCreate(PETSC_COMM_WORLD,&Kpf);CHKERRCONTINUE(ierr);
- ierr = MatCreate(PETSC_COMM_WORLD,&Kpp);CHKERRCONTINUE(ierr);
+// PetscErrorCode ierr = MatCreate(PETSC_COMM_WORLD,&Kpf);CHKERRCONTINUE(ierr);
+// ierr = MatCreate(PETSC_COMM_WORLD,&Kpp);CHKERRCONTINUE(ierr);
-#ifdef CPPUTILS_VERBOSE
- out<<" Allocating memory... "<<timer<<endl;
-#endif
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Allocating memory... "<<timer<<endl;
+// #endif
- ierr = MatSetFromOptions(Kpf);CHKERRCONTINUE(ierr);
- ierr = MatSetFromOptions(Kpp);CHKERRCONTINUE(ierr);
+// ierr = MatSetFromOptions(Kpf);CHKERRCONTINUE(ierr);
+// ierr = MatSetFromOptions(Kpp);CHKERRCONTINUE(ierr);
- ierr = MatSetSizes(Kpf,PETSC_DECIDE,PETSC_DECIDE, KKpf.rows(), KKpf.columns());CHKERRCONTINUE(ierr);
- ierr = MatSetSizes(Kpp,PETSC_DECIDE,PETSC_DECIDE, KKpp.rows(), KKpp.columns());CHKERRCONTINUE(ierr);
+// ierr = MatSetSizes(Kpf,PETSC_DECIDE,PETSC_DECIDE, KKpf.rows(), KKpf.columns());CHKERRCONTINUE(ierr);
+// ierr = MatSetSizes(Kpp,PETSC_DECIDE,PETSC_DECIDE, KKpp.rows(), KKpp.columns());CHKERRCONTINUE(ierr);
- // get number of non-zeros in both diagonal and non-diagonal portions of the matrix
+// // get number of non-zeros in both diagonal and non-diagonal portions of the matrix
- std::pair<size_t,size_t> Kpf_nz = KKpf.non_zero_off_diagonal();
- std::pair<size_t,size_t> Kpp_nz = KKpp.non_zero_off_diagonal();
+// std::pair<size_t,size_t> Kpf_nz = KKpf.non_zero_off_diagonal();
+// std::pair<size_t,size_t> Kpp_nz = KKpp.non_zero_off_diagonal();
- ierr = MatMPIAIJSetPreallocation(Kpf, Kpf_nz.first, PETSC_NULL, Kpf_nz.second, PETSC_NULL); CHKERRCONTINUE(ierr);
- ierr = MatMPIAIJSetPreallocation(Kpp, Kpp_nz.first, PETSC_NULL, Kpp_nz.second, PETSC_NULL); CHKERRCONTINUE(ierr);
- ierr = MatSeqAIJSetPreallocation(Kpf, Kpf_nz.first, PETSC_NULL); CHKERRCONTINUE(ierr);
- ierr = MatSeqAIJSetPreallocation(Kpp, Kpp_nz.first, PETSC_NULL); CHKERRCONTINUE(ierr);
+// ierr = MatMPIAIJSetPreallocation(Kpf, Kpf_nz.first, PETSC_NULL, Kpf_nz.second, PETSC_NULL); CHKERRCONTINUE(ierr);
+// ierr = MatMPIAIJSetPreallocation(Kpp, Kpp_nz.first, PETSC_NULL, Kpp_nz.second, PETSC_NULL); CHKERRCONTINUE(ierr);
+// ierr = MatSeqAIJSetPreallocation(Kpf, Kpf_nz.first, PETSC_NULL); CHKERRCONTINUE(ierr);
+// ierr = MatSeqAIJSetPreallocation(Kpp, Kpp_nz.first, PETSC_NULL); CHKERRCONTINUE(ierr);
- for (sparse_matrix_type::const_hash_iterator it = KKpf.map_.begin(); it != KKpf.map_.end(); ++it) {
+// for (sparse_matrix_type::const_hash_iterator it = KKpf.map_.begin(); it != KKpf.map_.end(); ++it) {
- // get subscripts
- std::pair<size_t,size_t> subs = KKpf.unhash(it->first);
- int row = subs.first;
- int col = subs.second;
- ierr = MatSetValues(Kpf, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERRCONTINUE(ierr);
- }
-
- for (sparse_matrix_type::const_hash_iterator it = KKpp.map_.begin(); it != KKpp.map_.end(); ++it) {
+// // get subscripts
+// std::pair<size_t,size_t> subs = KKpf.unhash(it->first);
+// int row = subs.first;
+// int col = subs.second;
+// ierr = MatSetValues(Kpf, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERRCONTINUE(ierr);
+// }
+
+// for (sparse_matrix_type::const_hash_iterator it = KKpp.map_.begin(); it != KKpp.map_.end(); ++it) {
- // get subscripts
- std::pair<size_t,size_t> subs = KKpp.unhash(it->first);
- int row = subs.first;
- int col = subs.second;
- ierr = MatSetValues(Kpf, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERRCONTINUE(ierr);
- }
+// // get subscripts
+// std::pair<size_t,size_t> subs = KKpp.unhash(it->first);
+// int row = subs.first;
+// int col = subs.second;
+// ierr = MatSetValues(Kpf, 1, &row, 1, &col, &it->second, ADD_VALUES);CHKERRCONTINUE(ierr);
+// }
-#ifdef CPPUTILS_VERBOSE
- out<<" Filled sparse matrices... "<<timer<<endl;
-#endif
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Filled sparse matrices... "<<timer<<endl;
+// #endif
- /*
- Assemble matrix, using the 2-step process:
- MatAssemblyBegin(), MatAssemblyEnd()
- Computations can be done while messages are in transition
- by placing code between these two statements.
- */
- ierr = MatAssemblyBegin(Kpf,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
- ierr = MatAssemblyBegin(Kpp,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
- ierr = MatAssemblyEnd(Kpf,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
- ierr = MatAssemblyEnd(Kpp,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
+// /*
+// Assemble matrix, using the 2-step process:
+// MatAssemblyBegin(), MatAssemblyEnd()
+// Computations can be done while messages are in transition
+// by placing code between these two statements.
+// */
+// ierr = MatAssemblyBegin(Kpf,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
+// ierr = MatAssemblyBegin(Kpp,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
+// ierr = MatAssemblyEnd(Kpf,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
+// ierr = MatAssemblyEnd(Kpp,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
-#ifdef CPPUTILS_VERBOSE
- out<<" Sparse matrices assembled... "<<timer<<endl;
-#endif
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Sparse matrices assembled... "<<timer<<endl;
+// #endif
- ierr = VecCreate(PETSC_COMM_WORLD,&Up);CHKERRCONTINUE(ierr);
- ierr = VecSetSizes(Up,PETSC_DECIDE, UUp.size());CHKERRCONTINUE(ierr);
- ierr = VecSetFromOptions(Up);CHKERRCONTINUE(ierr);
+// ierr = VecCreate(PETSC_COMM_WORLD,&Up);CHKERRCONTINUE(ierr);
+// ierr = VecSetSizes(Up,PETSC_DECIDE, UUp.size());CHKERRCONTINUE(ierr);
+// ierr = VecSetFromOptions(Up);CHKERRCONTINUE(ierr);
- ierr = VecCreate(PETSC_COMM_WORLD,&Pf);CHKERRCONTINUE(ierr);
- ierr = VecSetSizes(Pf,PETSC_DECIDE, KKpf.rows());CHKERRCONTINUE(ierr);
- ierr = VecSetFromOptions(Pf);CHKERRCONTINUE(ierr);
- ierr = VecDuplicate(Pf,&Pp);CHKERRCONTINUE(ierr);
+// ierr = VecCreate(PETSC_COMM_WORLD,&Pf);CHKERRCONTINUE(ierr);
+// ierr = VecSetSizes(Pf,PETSC_DECIDE, KKpf.rows());CHKERRCONTINUE(ierr);
+// ierr = VecSetFromOptions(Pf);CHKERRCONTINUE(ierr);
+// ierr = VecDuplicate(Pf,&Pp);CHKERRCONTINUE(ierr);
-#ifdef CPPUTILS_VERBOSE
- out<<" Vectors created... "<<timer<<endl;
-#endif
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Vectors created... "<<timer<<endl;
+// #endif
- /* Set hight-hand-side vector */
- for (sparse_vector_type::const_hash_iterator it = UUp.map_.begin(); it != UUp.map_.end(); ++it) {
- int row = it->first;
- ierr = VecSetValues(Up, 1, &row, &it->second, ADD_VALUES);
- }
+// /* Set hight-hand-side vector */
+// for (sparse_vector_type::const_hash_iterator it = UUp.map_.begin(); it != UUp.map_.end(); ++it) {
+// int row = it->first;
+// ierr = VecSetValues(Up, 1, &row, &it->second, ADD_VALUES);
+// }
- /*
- Assemble vector, using the 2-step process:
- VecAssemblyBegin(), VecAssemblyEnd()
- Computations can be done while messages are in transition
- by placing code between these two statements.
- */
- ierr = VecAssemblyBegin(Up);CHKERRCONTINUE(ierr);
- ierr = VecAssemblyEnd(Up);CHKERRCONTINUE(ierr);
+// /*
+// Assemble vector, using the 2-step process:
+// VecAssemblyBegin(), VecAssemblyEnd()
+// Computations can be done while messages are in transition
+// by placing code between these two statements.
+// */
+// ierr = VecAssemblyBegin(Up);CHKERRCONTINUE(ierr);
+// ierr = VecAssemblyEnd(Up);CHKERRCONTINUE(ierr);
- // add Kpf*Uf
- ierr = MatMult(Kpf, x_, Pf);
+// // add Kpf*Uf
+// ierr = MatMult(Kpf, x_, Pf);
- // add Kpp*Up
- ierr = MatMultAdd(Kpp, Up, Pf, Pp);
+// // add Kpp*Up
+// ierr = MatMultAdd(Kpp, Up, Pf, Pp);
-#ifdef CPPUTILS_VERBOSE
- out<<" Matrices multiplied... "<<timer<<endl;
-#endif
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Matrices multiplied... "<<timer<<endl;
+// #endif
- VecScatter ctx;
- Vec Pp_all;
+// VecScatter ctx;
+// Vec Pp_all;
- ierr = VecScatterCreateToAll(Pp, &ctx, &Pp_all);CHKERRCONTINUE(ierr);
- ierr = VecScatterBegin(ctx,Pp,Pp_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
- ierr = VecScatterEnd(ctx,Pp,Pp_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
+// ierr = VecScatterCreateToAll(Pp, &ctx, &Pp_all);CHKERRCONTINUE(ierr);
+// ierr = VecScatterBegin(ctx,Pp,Pp_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
+// ierr = VecScatterEnd(ctx,Pp,Pp_all,INSERT_VALUES,SCATTER_FORWARD);CHKERRCONTINUE(ierr);
-#ifdef CPPUTILS_VERBOSE
- out<<" Vector scattered... "<<timer<<endl;
-#endif
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Vector scattered... "<<timer<<endl;
+// #endif
- /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- Get values from solution and store them in the object that will be
- returned
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+// /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+// Get values from solution and store them in the object that will be
+// returned
+// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
- sparse_vector_type pp(KKpf.rows());
+// sparse_vector_type pp(KKpf.rows());
- // get reaction vector
- for (int i=0; i<KKpf.rows(); ++i) {
+// // get reaction vector
+// for (int i=0; i<KKpf.rows(); ++i) {
- PetscScalar v;
- ierr = VecGetValues(Pp_all, 1, &i, &v);
- if (v != PetscScalar())
- pp[i] = v;
- }
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Vector copied... "<<timer<<endl;
-#endif
-
- ierr = MatDestroy(&Kpf);CHKERRCONTINUE(ierr);
- ierr = MatDestroy(&Kpp);CHKERRCONTINUE(ierr);
- ierr = VecDestroy(&Up);CHKERRCONTINUE(ierr);
- ierr = VecDestroy(&Pf);CHKERRCONTINUE(ierr);
- ierr = VecDestroy(&Pp);CHKERRCONTINUE(ierr);
- ierr = VecDestroy(&Pp_all);CHKERRCONTINUE(ierr);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Temporary data structures destroyed... "<<timer<<endl;
-#endif
-
- return pp;
-}
+// PetscScalar v;
+// ierr = VecGetValues(Pp_all, 1, &i, &v);
+// if (v != PetscScalar())
+// pp[i] = v;
+// }
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Vector copied... "<<timer<<endl;
+// #endif
+
+// ierr = MatDestroy(&Kpf);CHKERRCONTINUE(ierr);
+// ierr = MatDestroy(&Kpp);CHKERRCONTINUE(ierr);
+// ierr = VecDestroy(&Up);CHKERRCONTINUE(ierr);
+// ierr = VecDestroy(&Pf);CHKERRCONTINUE(ierr);
+// ierr = VecDestroy(&Pp);CHKERRCONTINUE(ierr);
+// ierr = VecDestroy(&Pp_all);CHKERRCONTINUE(ierr);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Temporary data structures destroyed... "<<timer<<endl;
+// #endif
+
+// return pp;
+// }
-SolverPETSc::sparse_vector_type SolverPETSc::operator()(const SolverPETSc::sparse_vector_type& aa, const SolverPETSc::sparse_vector_type& bb) {
-
- assert(aa.size() == bb.size());
-
-#ifdef CPPUTILS_VERBOSE
- // parallel output stream
- Output_stream out;
- // timer
- cpputils::ctimer timer;
- out<<"Inside SolverPETSc::operator()(const sparse_vector&, const sparse_vector&). "<<timer<<endl;
-#endif
-
- Vec r;
-
- PetscErrorCode ierr = VecCreate(PETSC_COMM_WORLD,&r);CHKERRCONTINUE(ierr);
- ierr = VecSetSizes(r,PETSC_DECIDE, aa.size());CHKERRCONTINUE(ierr);
- ierr = VecSetFromOptions(r);CHKERRCONTINUE(ierr);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Vectors created... "<<timer<<endl;
-#endif
-
- // set values
- for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != aa.map_.end(); ++it) {
- int row = it->first;
- ierr = VecSetValues(r, 1, &row, &it->second, ADD_VALUES);
- }
- for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
- int row = it->first;
- ierr = VecSetValues(r, 1, &row, &it->second, ADD_VALUES);
- }
-
- /*
- Assemble vector, using the 2-step process:
- VecAssemblyBegin(), VecAssemblyEnd()
- Computations can be done while messages are in transition
- by placing code between these two statements.
- */
- ierr = VecAssemblyBegin(r);CHKERRCONTINUE(ierr);
- ierr = VecAssemblyEnd(r);CHKERRCONTINUE(ierr);
-
-
- sparse_vector_type rr(aa.size());
-
- for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != aa.map_.end(); ++it) {
- int row = it->first;
- ierr = VecGetValues(r, 1, &row, &rr[row]);
- }
- for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
- int row = it->first;
- ierr = VecGetValues(r, 1, &row, &rr[row]);
- }
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Vector copied... "<<timer<<endl;
-#endif
-
- ierr = VecDestroy(&r);CHKERRCONTINUE(ierr);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Temporary data structures destroyed... "<<timer<<endl;
-#endif
-
- return rr;
-}
+// SolverPETSc::sparse_vector_type SolverPETSc::operator()(const SolverPETSc::sparse_vector_type& aa, const SolverPETSc::sparse_vector_type& bb) {
+
+// assert(aa.size() == bb.size());
+
+// #ifdef CPPUTILS_VERBOSE
+// // parallel output stream
+// Output_stream out;
+// // timer
+// cpputils::ctimer timer;
+// out<<"Inside SolverPETSc::operator()(const sparse_vector&, const sparse_vector&). "<<timer<<endl;
+// #endif
+
+// Vec r;
+
+// PetscErrorCode ierr = VecCreate(PETSC_COMM_WORLD,&r);CHKERRCONTINUE(ierr);
+// ierr = VecSetSizes(r,PETSC_DECIDE, aa.size());CHKERRCONTINUE(ierr);
+// ierr = VecSetFromOptions(r);CHKERRCONTINUE(ierr);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Vectors created... "<<timer<<endl;
+// #endif
+
+// // set values
+// for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != aa.map_.end(); ++it) {
+// int row = it->first;
+// ierr = VecSetValues(r, 1, &row, &it->second, ADD_VALUES);
+// }
+// for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
+// int row = it->first;
+// ierr = VecSetValues(r, 1, &row, &it->second, ADD_VALUES);
+// }
+
+// /*
+// Assemble vector, using the 2-step process:
+// VecAssemblyBegin(), VecAssemblyEnd()
+// Computations can be done while messages are in transition
+// by placing code between these two statements.
+// */
+// ierr = VecAssemblyBegin(r);CHKERRCONTINUE(ierr);
+// ierr = VecAssemblyEnd(r);CHKERRCONTINUE(ierr);
+
+
+// sparse_vector_type rr(aa.size());
+
+// for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != aa.map_.end(); ++it) {
+// int row = it->first;
+// ierr = VecGetValues(r, 1, &row, &rr[row]);
+// }
+// for (sparse_vector_type::const_hash_iterator it = bb.map_.begin(); it != bb.map_.end(); ++it) {
+// int row = it->first;
+// ierr = VecGetValues(r, 1, &row, &rr[row]);
+// }
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Vector copied... "<<timer<<endl;
+// #endif
+
+// ierr = VecDestroy(&r);CHKERRCONTINUE(ierr);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Temporary data structures destroyed... "<<timer<<endl;
+// #endif
+
+// return rr;
+// }
-SolverPETSc::value_type SolverPETSc::norm(const SolverPETSc::sparse_matrix_type& aa, Element_insertion_type flag) {
-
-#ifdef CPPUTILS_VERBOSE
- // parallel output stream
- Output_stream out;
- // timer
- cpputils::ctimer timer;
- out<<"Inside SolverPETSc::norm(const sparse_matrix&). "<<timer<<endl;
-#endif
-
- Mat r;
-
- PetscErrorCode ierr = MatCreate(PETSC_COMM_WORLD,&r);CHKERRCONTINUE(ierr);
- ierr = MatSetSizes(r,PETSC_DECIDE,PETSC_DECIDE, aa.rows(), aa.columns());CHKERRCONTINUE(ierr);
- ierr = MatSetFromOptions(r);CHKERRCONTINUE(ierr);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Matrix created... "<<timer<<endl;
-#endif
-
- // set values
- for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != aa.map_.end(); ++it) {
- // get subscripts
- std::pair<size_t,size_t> subs = aa.unhash(it->first);
- int row = subs.first;
- int col = subs.second;
+// SolverPETSc::value_type SolverPETSc::norm(const SolverPETSc::sparse_matrix_type& aa, Element_insertion_type flag) {
+
+// #ifdef CPPUTILS_VERBOSE
+// // parallel output stream
+// Output_stream out;
+// // timer
+// cpputils::ctimer timer;
+// out<<"Inside SolverPETSc::norm(const sparse_matrix&). "<<timer<<endl;
+// #endif
+
+// Mat r;
+
+// PetscErrorCode ierr = MatCreate(PETSC_COMM_WORLD,&r);CHKERRCONTINUE(ierr);
+// ierr = MatSetSizes(r,PETSC_DECIDE,PETSC_DECIDE, aa.rows(), aa.columns());CHKERRCONTINUE(ierr);
+// ierr = MatSetFromOptions(r);CHKERRCONTINUE(ierr);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Matrix created... "<<timer<<endl;
+// #endif
+
+// // set values
+// for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != aa.map_.end(); ++it) {
+// // get subscripts
+// std::pair<size_t,size_t> subs = aa.unhash(it->first);
+// int row = subs.first;
+// int col = subs.second;
- if (flag == Add_t)
- ierr = MatSetValues(r, 1, &row, 1, &col, &it->second, ADD_VALUES);
- else if (flag == Insert_t)
- ierr = MatSetValues(r, 1, &row, 1, &col, &it->second, INSERT_VALUES);
- CHKERRCONTINUE(ierr);
- }
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Matrix filled..."<<timer<<endl;
-#endif
-
- /*
- Assemble vector, using the 2-step process:
- VecAssemblyBegin(), VecAssemblyEnd()
- Computations can be done while messages are in transition
- by placing code between these two statements.
- */
- ierr = MatAssemblyBegin(r,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
- ierr = MatAssemblyEnd(r,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
-
- value_type nrm;
-
- MatNorm(r,NORM_FROBENIUS,&nrm);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Norm computed... "<<timer<<endl;
-#endif
-
- ierr = MatDestroy(&r);CHKERRCONTINUE(ierr);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Temporary data structures destroyed... "<<timer<<endl;
-#endif
-
- return nrm;
-}
+// if (flag == Add_t)
+// ierr = MatSetValues(r, 1, &row, 1, &col, &it->second, ADD_VALUES);
+// else if (flag == Insert_t)
+// ierr = MatSetValues(r, 1, &row, 1, &col, &it->second, INSERT_VALUES);
+// CHKERRCONTINUE(ierr);
+// }
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Matrix filled..."<<timer<<endl;
+// #endif
+
+// /*
+// Assemble vector, using the 2-step process:
+// VecAssemblyBegin(), VecAssemblyEnd()
+// Computations can be done while messages are in transition
+// by placing code between these two statements.
+// */
+// ierr = MatAssemblyBegin(r,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
+// ierr = MatAssemblyEnd(r,MAT_FINAL_ASSEMBLY);CHKERRCONTINUE(ierr);
+
+// value_type nrm;
+
+// MatNorm(r,NORM_FROBENIUS,&nrm);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Norm computed... "<<timer<<endl;
+// #endif
+
+// ierr = MatDestroy(&r);CHKERRCONTINUE(ierr);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Temporary data structures destroyed... "<<timer<<endl;
+// #endif
+
+// return nrm;
+// }
-SolverPETSc::value_type SolverPETSc::norm(const SolverPETSc::sparse_vector_type& aa, Element_insertion_type flag) {
-
-#ifdef CPPUTILS_VERBOSE
- // parallel output stream
- Output_stream out;
- // timer
- cpputils::ctimer timer;
- out<<"Inside SolverPETSc::norm(const sparse_vector&). "<<timer<<endl;
-#endif
-
- Vec r;
-
- PetscErrorCode ierr = VecCreate(PETSC_COMM_WORLD,&r);CHKERRCONTINUE(ierr);
- ierr = VecSetSizes(r,PETSC_DECIDE, aa.size());CHKERRCONTINUE(ierr);
- ierr = VecSetFromOptions(r);CHKERRCONTINUE(ierr);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Vector created... "<<timer<<endl;
-#endif
-
- // set values
- for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != aa.map_.end(); ++it) {
- int row = it->first;
- if (flag == Add_t)
- ierr = VecSetValues(r, 1, &row, &it->second, ADD_VALUES);
- else if (flag == Insert_t)
- ierr = VecSetValues(r, 1, &row, &it->second, INSERT_VALUES);
- CHKERRCONTINUE(ierr);
- }
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Vector filled..."<<timer<<endl;
-#endif
-
- /*
- Assemble vector, using the 2-step process:
- VecAssemblyBegin(), VecAssemblyEnd()
- Computations can be done while messages are in transition
- by placing code between these two statements.
- */
- ierr = VecAssemblyBegin(r);CHKERRCONTINUE(ierr);
- ierr = VecAssemblyEnd(r);CHKERRCONTINUE(ierr);
-
- value_type nrm;
-
- VecNorm(r,NORM_2,&nrm);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Norm computed... "<<timer<<endl;
-#endif
-
- ierr = VecDestroy(&r);CHKERRCONTINUE(ierr);
-
-#ifdef CPPUTILS_VERBOSE
- out<<" Temporary data structures destroyed... "<<timer<<endl;
-#endif
-
- return nrm;
-
-}
+// SolverPETSc::value_type SolverPETSc::norm(const SolverPETSc::sparse_vector_type& aa, Element_insertion_type flag) {
+
+// #ifdef CPPUTILS_VERBOSE
+// // parallel output stream
+// Output_stream out;
+// // timer
+// cpputils::ctimer timer;
+// out<<"Inside SolverPETSc::norm(const sparse_vector&). "<<timer<<endl;
+// #endif
+
+// Vec r;
+
+// PetscErrorCode ierr = VecCreate(PETSC_COMM_WORLD,&r);CHKERRCONTINUE(ierr);
+// ierr = VecSetSizes(r,PETSC_DECIDE, aa.size());CHKERRCONTINUE(ierr);
+// ierr = VecSetFromOptions(r);CHKERRCONTINUE(ierr);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Vector created... "<<timer<<endl;
+// #endif
+
+// // set values
+// for (sparse_vector_type::const_hash_iterator it = aa.map_.begin(); it != aa.map_.end(); ++it) {
+// int row = it->first;
+// if (flag == Add_t)
+// ierr = VecSetValues(r, 1, &row, &it->second, ADD_VALUES);
+// else if (flag == Insert_t)
+// ierr = VecSetValues(r, 1, &row, &it->second, INSERT_VALUES);
+// CHKERRCONTINUE(ierr);
+// }
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Vector filled..."<<timer<<endl;
+// #endif
+
+// /*
+// Assemble vector, using the 2-step process:
+// VecAssemblyBegin(), VecAssemblyEnd()
+// Computations can be done while messages are in transition
+// by placing code between these two statements.
+// */
+// ierr = VecAssemblyBegin(r);CHKERRCONTINUE(ierr);
+// ierr = VecAssemblyEnd(r);CHKERRCONTINUE(ierr);
+
+// value_type nrm;
+
+// VecNorm(r,NORM_2,&nrm);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Norm computed... "<<timer<<endl;
+// #endif
+
+// ierr = VecDestroy(&r);CHKERRCONTINUE(ierr);
+
+// #ifdef CPPUTILS_VERBOSE
+// out<<" Temporary data structures destroyed... "<<timer<<endl;
+// #endif
+
+// return nrm;
+
+// }
-//
-///* -------------------------------------------------------------------------- */
-//SolverMumps::SolverMumps(SparseMatrix & matrix,
-// const ID & id,
-// const MemoryID & memory_id) :
-//Solver(matrix, id, memory_id), is_mumps_data_initialized(false), rhs_is_local(true) {
-// AKANTU_DEBUG_IN();
-//
-//#ifdef AKANTU_USE_MPI
-// parallel_method = SolverMumpsOptions::_fully_distributed;
-//#else //AKANTU_USE_MPI
-// parallel_method = SolverMumpsOptions::_master_slave_distributed;
-//#endif //AKANTU_USE_MPI
-//
-// CommunicatorEventHandler & comm_event_handler = *this;
-//
-// communicator.registerEventHandler(comm_event_handler);
-//
-// AKANTU_DEBUG_OUT();
-//}
-//
-///* -------------------------------------------------------------------------- */
-//SolverMumps::~SolverMumps() {
-// AKANTU_DEBUG_IN();
-//
-// AKANTU_DEBUG_OUT();
-//}
-//
-///* -------------------------------------------------------------------------- */
-//void SolverMumps::destroyMumpsData() {
-// AKANTU_DEBUG_IN();
-//
-// if(is_mumps_data_initialized) {
-// mumps_data.job = _smj_destroy; // destroy
-// dmumps_c(&mumps_data);
-// is_mumps_data_initialized = false;
-// }
-//
-// AKANTU_DEBUG_OUT();
-//}
-//
-///* -------------------------------------------------------------------------- */
-//void SolverMumps::onCommunicatorFinalize(const StaticCommunicator & comm) {
-// AKANTU_DEBUG_IN();
-//
-// try{
-// const StaticCommunicatorMPI & comm_mpi =
-// dynamic_cast<const StaticCommunicatorMPI &>(comm.getRealStaticCommunicator());
-// if(mumps_data.comm_fortran == MPI_Comm_c2f(comm_mpi.getMPICommunicator()))
-// destroyMumpsData();
-// } catch(...) {}
-//
-// AKANTU_DEBUG_OUT();
-//}
-//
-///* -------------------------------------------------------------------------- */
-//void SolverMumps::initMumpsData(SolverMumpsOptions::ParallelMethod parallel_method) {
-// switch(parallel_method) {
-// case SolverMumpsOptions::_fully_distributed:
-// icntl(18) = 3; //fully distributed
-// icntl(28) = 0; //automatic choice
-//
-// mumps_data.nz_loc = matrix->getNbNonZero();
-// mumps_data.irn_loc = matrix->getIRN().values;
-// mumps_data.jcn_loc = matrix->getJCN().values;
-// break;
-// case SolverMumpsOptions::_master_slave_distributed:
-// if(prank == 0) {
-// mumps_data.nz = matrix->getNbNonZero();
-// mumps_data.irn = matrix->getIRN().values;
-// mumps_data.jcn = matrix->getJCN().values;
-// } else {
-// mumps_data.nz = 0;
-// mumps_data.irn = NULL;
-// mumps_data.jcn = NULL;
-//
-// icntl(18) = 0; //centralized
-// icntl(28) = 0; //sequential analysis
-// }
-// break;
-// }
-//}
-//
-///* -------------------------------------------------------------------------- */
-//void SolverMumps::initialize(SolverOptions & options) {
-// AKANTU_DEBUG_IN();
-//
-// mumps_data.par = 1;
-//
-// if(SolverMumpsOptions * opt = dynamic_cast<SolverMumpsOptions *>(&options)) {
-// if(opt->parallel_method == SolverMumpsOptions::_master_slave_distributed) {
-// mumps_data.par = 0;
-// }
-// }
-//
-// mumps_data.sym = 2 * (matrix->getSparseMatrixType() == _symmetric);
-// prank = communicator.whoAmI();
-//#ifdef AKANTU_USE_MPI
-// mumps_data.comm_fortran = MPI_Comm_c2f(dynamic_cast<const StaticCommunicatorMPI &>(communicator.getRealStaticCommunicator()).getMPICommunicator());
-//#endif
-//
-// if(AKANTU_DEBUG_TEST(dblTrace)) {
-// icntl(1) = 2;
-// icntl(2) = 2;
-// icntl(3) = 2;
-// icntl(4) = 4;
-// }
-//
-// mumps_data.job = _smj_initialize; //initialize
-// dmumps_c(&mumps_data);
-// is_mumps_data_initialized = true;
-//
-// /* ------------------------------------------------------------------------ */
-// UInt size = matrix->getSize();
-//
-// if(prank == 0) {
-// std::stringstream sstr_rhs; sstr_rhs << id << ":rhs";
-// rhs = &(alloc<Real>(sstr_rhs.str(), size, 1, REAL_INIT_VALUE));
-// } else {
-// rhs = NULL;
-// }
-//
-// /// No outputs
-// icntl(1) = 0;
-// icntl(2) = 0;
-// icntl(3) = 0;
-// icntl(4) = 0;
-// mumps_data.nz_alloc = 0;
-//
-// if (AKANTU_DEBUG_TEST(dblDump)) icntl(4) = 4;
-//
-// mumps_data.n = size;
-//
-// if(AKANTU_DEBUG_TEST(dblDump)) {
-// strcpy(mumps_data.write_problem, "mumps_matrix.mtx");
-// }
-//
-// /* ------------------------------------------------------------------------ */
-// // Default Scaling
-// icntl(8) = 77;
-//
-// icntl(5) = 0; // Assembled matrix
-//
-// SolverMumpsOptions * opt = dynamic_cast<SolverMumpsOptions *>(&options);
-// if(opt)
-// parallel_method = opt->parallel_method;
-//
-// initMumpsData(parallel_method);
-//
-// mumps_data.job = _smj_analyze; //analyze
-// dmumps_c(&mumps_data);
-//
-// AKANTU_DEBUG_OUT();
-//}
-//
-///* -------------------------------------------------------------------------- */
-//void SolverMumps::setRHS(Array<Real> & rhs) {
-// if(prank == 0) {
-// matrix->getDOFSynchronizer().gather(rhs, 0, this->rhs);
-// } else {
-// matrix->getDOFSynchronizer().gather(rhs, 0);
-// }
-//}
-//
-///* -------------------------------------------------------------------------- */
-//void SolverMumps::solve() {
-// AKANTU_DEBUG_IN();
-//
-// if(parallel_method == SolverMumpsOptions::_fully_distributed)
-// mumps_data.a_loc = matrix->getA().values;
-// else
-// if(prank == 0) {
-// mumps_data.a = matrix->getA().values;
-// }
-//
-// if(prank == 0) {
-// mumps_data.rhs = rhs->values;
-// }
-//
-// /// Default centralized dense second member
-// icntl(20) = 0;
-// icntl(21) = 0;
-//
-// mumps_data.job = _smj_factorize_solve; //solve
-// dmumps_c(&mumps_data);
-//
-// AKANTU_DEBUG_ASSERT(info(1) != -10, "Singular matrix");
-// AKANTU_DEBUG_ASSERT(info(1) == 0,
-// "Error in mumps during solve process, check mumps user guide INFO(1) ="
-// << info(1));
-//
-// AKANTU_DEBUG_OUT();
-//}
-//
-///* -------------------------------------------------------------------------- */
-//void SolverMumps::solve(Array<Real> & solution) {
-// AKANTU_DEBUG_IN();
-//
-// solve();
-//
-// if(prank == 0) {
-// matrix->getDOFSynchronizer().scatter(solution, 0, this->rhs);
-// } else {
-// matrix->getDOFSynchronizer().scatter(solution, 0);
-// }
-//
-// AKANTU_DEBUG_OUT();
-//}
+// //
+// ///* -------------------------------------------------------------------------- */
+// //SolverMumps::SolverMumps(SparseMatrix & matrix,
+// // const ID & id,
+// // const MemoryID & memory_id) :
+// //Solver(matrix, id, memory_id), is_mumps_data_initialized(false), rhs_is_local(true) {
+// // AKANTU_DEBUG_IN();
+// //
+// //#ifdef AKANTU_USE_MPI
+// // parallel_method = SolverMumpsOptions::_fully_distributed;
+// //#else //AKANTU_USE_MPI
+// // parallel_method = SolverMumpsOptions::_master_slave_distributed;
+// //#endif //AKANTU_USE_MPI
+// //
+// // CommunicatorEventHandler & comm_event_handler = *this;
+// //
+// // communicator.registerEventHandler(comm_event_handler);
+// //
+// // AKANTU_DEBUG_OUT();
+// //}
+// //
+// ///* -------------------------------------------------------------------------- */
+// //SolverMumps::~SolverMumps() {
+// // AKANTU_DEBUG_IN();
+// //
+// // AKANTU_DEBUG_OUT();
+// //}
+// //
+// ///* -------------------------------------------------------------------------- */
+// //void SolverMumps::destroyMumpsData() {
+// // AKANTU_DEBUG_IN();
+// //
+// // if(is_mumps_data_initialized) {
+// // mumps_data.job = _smj_destroy; // destroy
+// // dmumps_c(&mumps_data);
+// // is_mumps_data_initialized = false;
+// // }
+// //
+// // AKANTU_DEBUG_OUT();
+// //}
+// //
+// ///* -------------------------------------------------------------------------- */
+// //void SolverMumps::onCommunicatorFinalize(const StaticCommunicator & comm) {
+// // AKANTU_DEBUG_IN();
+// //
+// // try{
+// // const StaticCommunicatorMPI & comm_mpi =
+// // dynamic_cast<const StaticCommunicatorMPI &>(comm.getRealStaticCommunicator());
+// // if(mumps_data.comm_fortran == MPI_Comm_c2f(comm_mpi.getMPICommunicator()))
+// // destroyMumpsData();
+// // } catch(...) {}
+// //
+// // AKANTU_DEBUG_OUT();
+// //}
+// //
+// ///* -------------------------------------------------------------------------- */
+// //void SolverMumps::initMumpsData(SolverMumpsOptions::ParallelMethod parallel_method) {
+// // switch(parallel_method) {
+// // case SolverMumpsOptions::_fully_distributed:
+// // icntl(18) = 3; //fully distributed
+// // icntl(28) = 0; //automatic choice
+// //
+// // mumps_data.nz_loc = matrix->getNbNonZero();
+// // mumps_data.irn_loc = matrix->getIRN().values;
+// // mumps_data.jcn_loc = matrix->getJCN().values;
+// // break;
+// // case SolverMumpsOptions::_master_slave_distributed:
+// // if(prank == 0) {
+// // mumps_data.nz = matrix->getNbNonZero();
+// // mumps_data.irn = matrix->getIRN().values;
+// // mumps_data.jcn = matrix->getJCN().values;
+// // } else {
+// // mumps_data.nz = 0;
+// // mumps_data.irn = NULL;
+// // mumps_data.jcn = NULL;
+// //
+// // icntl(18) = 0; //centralized
+// // icntl(28) = 0; //sequential analysis
+// // }
+// // break;
+// // }
+// //}
+// //
+// ///* -------------------------------------------------------------------------- */
+// //void SolverMumps::initialize(SolverOptions & options) {
+// // AKANTU_DEBUG_IN();
+// //
+// // mumps_data.par = 1;
+// //
+// // if(SolverMumpsOptions * opt = dynamic_cast<SolverMumpsOptions *>(&options)) {
+// // if(opt->parallel_method == SolverMumpsOptions::_master_slave_distributed) {
+// // mumps_data.par = 0;
+// // }
+// // }
+// //
+// // mumps_data.sym = 2 * (matrix->getSparseMatrixType() == _symmetric);
+// // prank = communicator.whoAmI();
+// //#ifdef AKANTU_USE_MPI
+// // mumps_data.comm_fortran = MPI_Comm_c2f(dynamic_cast<const StaticCommunicatorMPI &>(communicator.getRealStaticCommunicator()).getMPICommunicator());
+// //#endif
+// //
+// // if(AKANTU_DEBUG_TEST(dblTrace)) {
+// // icntl(1) = 2;
+// // icntl(2) = 2;
+// // icntl(3) = 2;
+// // icntl(4) = 4;
+// // }
+// //
+// // mumps_data.job = _smj_initialize; //initialize
+// // dmumps_c(&mumps_data);
+// // is_mumps_data_initialized = true;
+// //
+// // /* ------------------------------------------------------------------------ */
+// // UInt size = matrix->getSize();
+// //
+// // if(prank == 0) {
+// // std::stringstream sstr_rhs; sstr_rhs << id << ":rhs";
+// // rhs = &(alloc<Real>(sstr_rhs.str(), size, 1, REAL_INIT_VALUE));
+// // } else {
+// // rhs = NULL;
+// // }
+// //
+// // /// No outputs
+// // icntl(1) = 0;
+// // icntl(2) = 0;
+// // icntl(3) = 0;
+// // icntl(4) = 0;
+// // mumps_data.nz_alloc = 0;
+// //
+// // if (AKANTU_DEBUG_TEST(dblDump)) icntl(4) = 4;
+// //
+// // mumps_data.n = size;
+// //
+// // if(AKANTU_DEBUG_TEST(dblDump)) {
+// // strcpy(mumps_data.write_problem, "mumps_matrix.mtx");
+// // }
+// //
+// // /* ------------------------------------------------------------------------ */
+// // // Default Scaling
+// // icntl(8) = 77;
+// //
+// // icntl(5) = 0; // Assembled matrix
+// //
+// // SolverMumpsOptions * opt = dynamic_cast<SolverMumpsOptions *>(&options);
+// // if(opt)
+// // parallel_method = opt->parallel_method;
+// //
+// // initMumpsData(parallel_method);
+// //
+// // mumps_data.job = _smj_analyze; //analyze
+// // dmumps_c(&mumps_data);
+// //
+// // AKANTU_DEBUG_OUT();
+// //}
+// //
+// ///* -------------------------------------------------------------------------- */
+// //void SolverMumps::setRHS(Array<Real> & rhs) {
+// // if(prank == 0) {
+// // matrix->getDOFSynchronizer().gather(rhs, 0, this->rhs);
+// // } else {
+// // matrix->getDOFSynchronizer().gather(rhs, 0);
+// // }
+// //}
+// //
+// ///* -------------------------------------------------------------------------- */
+// //void SolverMumps::solve() {
+// // AKANTU_DEBUG_IN();
+// //
+// // if(parallel_method == SolverMumpsOptions::_fully_distributed)
+// // mumps_data.a_loc = matrix->getA().values;
+// // else
+// // if(prank == 0) {
+// // mumps_data.a = matrix->getA().values;
+// // }
+// //
+// // if(prank == 0) {
+// // mumps_data.rhs = rhs->values;
+// // }
+// //
+// // /// Default centralized dense second member
+// // icntl(20) = 0;
+// // icntl(21) = 0;
+// //
+// // mumps_data.job = _smj_factorize_solve; //solve
+// // dmumps_c(&mumps_data);
+// //
+// // AKANTU_DEBUG_ASSERT(info(1) != -10, "Singular matrix");
+// // AKANTU_DEBUG_ASSERT(info(1) == 0,
+// // "Error in mumps during solve process, check mumps user guide INFO(1) ="
+// // << info(1));
+// //
+// // AKANTU_DEBUG_OUT();
+// //}
+// //
+// ///* -------------------------------------------------------------------------- */
+// //void SolverMumps::solve(Array<Real> & solution) {
+// // AKANTU_DEBUG_IN();
+// //
+// // solve();
+// //
+// // if(prank == 0) {
+// // matrix->getDOFSynchronizer().scatter(solution, 0, this->rhs);
+// // } else {
+// // matrix->getDOFSynchronizer().scatter(solution, 0);
+// // }
+// //
+// // AKANTU_DEBUG_OUT();
+// //}
-__END_AKANTU__
+// __END_AKANTU__
diff --git a/src/solver/solver_petsc.hh b/src/solver/solver_petsc.hh
index f4b7e1c42..352a2db88 100644
--- a/src/solver/solver_petsc.hh
+++ b/src/solver/solver_petsc.hh
@@ -1,273 +1,271 @@
/**
* @file solver_petsc.hh
*
# @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Mon Dec 13 10:48:06 2010
*
* @brief Solver class implementation for the petsc solver
*
* @section LICENSE
*
* Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __AKANTU_SOLVER_PETSC_HH__
-#define __AKANTU_SOLVER_PETSC_HH__
+// #ifndef __AKANTU_SOLVER_PETSC_HH__
+// #define __AKANTU_SOLVER_PETSC_HH__
-#include <petscksp.h>
+// #include <petscksp.h>
-#include "solver.hh"
-#include "static_communicator.hh"
-#include "sparse_matrix.hh"
+// #include "solver.hh"
+// #include "static_communicator.hh"
+// #include "sparse_matrix.hh"
-struct Mat;
+// __BEGIN_AKANTU__
-__BEGIN_AKANTU__
-
-struct SolverPETSc : public Solver, public CommunicatorEventHandler {
+// struct SolverPETSc : public Solver, public CommunicatorEventHandler {
- typedef double value_type;
- typedef sparse_vector<value_type> sparse_vector_type;
- typedef SparseMatrix sparse_matrix_type;
+// typedef double value_type;
+// typedef sparse_vector<value_type> sparse_vector_type;
+// typedef SparseMatrix sparse_matrix_type;
- Mat A_; //!< linear system matrix
- Vec x_; //!< Solution vector
- KSP ksp_; //!< linear solver context
+// Mat A_; //!< linear system matrix
+// Vec x_; //!< Solution vector
+// KSP ksp_; //!< linear solver context
- bool allocated_;
+// bool allocated_;
- SolverPETSc(int argc, char *argv[]) : allocated_(false) {
+// SolverPETSc(int argc, char *argv[]) : allocated_(false) {
- PetscInitialize(&argc, &argv,NULL,NULL);
- PetscErrorCode ierr;
+// PetscInitialize(&argc, &argv,NULL,NULL);
+// PetscErrorCode ierr;
- // create linear solver context
- ierr = KSPCreate(PETSC_COMM_WORLD, &ksp_);CHKERRCONTINUE(ierr);
+// // create linear solver context
+// ierr = KSPCreate(PETSC_COMM_WORLD, &ksp_);CHKERRCONTINUE(ierr);
- // initial nonzero guess
- ierr = KSPSetInitialGuessNonzero(ksp_,PETSC_TRUE); CHKERRCONTINUE(ierr);
+// // initial nonzero guess
+// ierr = KSPSetInitialGuessNonzero(ksp_,PETSC_TRUE); CHKERRCONTINUE(ierr);
- // set runtime options
- ierr = KSPSetFromOptions(ksp_);CHKERRCONTINUE(ierr);
+// // set runtime options
+// ierr = KSPSetFromOptions(ksp_);CHKERRCONTINUE(ierr);
- /*
- Set linear solver defaults for this problem (optional).
- - By extracting the KSP and PC contexts from the KSP context,
- we can then directly call any KSP and PC routines to set
- various options.
- - The following four statements are optional; all of these
- parameters could alternatively be specified at runtime via
- KSPSetFromOptions();
- */
- // ierr = KSPGetPC(ksp_,&pc);CHKERRCONTINUE(ierr);
- // ierr = PCSetType(pc,PCILU);CHKERRCONTINUE(ierr);
- // ierr = PCSetType(pc,PCJACOBI);CHKERRCONTINUE(ierr);
- ierr = KSPSetTolerances(ksp_,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRCONTINUE(ierr);
- }
-
- //! Overload operator() to solve system of linear equations
- sparse_vector_type operator()(const sparse_matrix_type& AA, const sparse_vector_type& bb);
-
- //! Overload operator() to obtain reaction vector
- sparse_vector_type operator()(const sparse_matrix_type& Kpf, const sparse_matrix_type& Kpp, const sparse_vector_type& Up);
-
- //! Overload operator() to obtain the addition two vectors
- sparse_vector_type operator()(const sparse_vector_type& aa, const sparse_vector_type& bb);
-
- value_type norm(const sparse_matrix_type& aa, Element_insertion_type it = Add_t);
-
- value_type norm(const sparse_vector_type& aa, Element_insertion_type it = Add_t);
-
- // NOTE: the destructor will return an error if it is called after MPI_Finalize is
- // called because it uses collect communication to free-up allocated memory.
- ~SolverPETSc() {
+// /*
+// Set linear solver defaults for this problem (optional).
+// - By extracting the KSP and PC contexts from the KSP context,
+// we can then directly call any KSP and PC routines to set
+// various options.
+// - The following four statements are optional; all of these
+// parameters could alternatively be specified at runtime via
+// KSPSetFromOptions();
+// */
+// // ierr = KSPGetPC(ksp_,&pc);CHKERRCONTINUE(ierr);
+// // ierr = PCSetType(pc,PCILU);CHKERRCONTINUE(ierr);
+// // ierr = PCSetType(pc,PCJACOBI);CHKERRCONTINUE(ierr);
+// ierr = KSPSetTolerances(ksp_,1.e-5,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRCONTINUE(ierr);
+// }
+
+// //! Overload operator() to solve system of linear equations
+// sparse_vector_type operator()(const sparse_matrix_type& AA, const sparse_vector_type& bb);
+
+// //! Overload operator() to obtain reaction vector
+// sparse_vector_type operator()(const sparse_matrix_type& Kpf, const sparse_matrix_type& Kpp, const sparse_vector_type& Up);
+
+// //! Overload operator() to obtain the addition two vectors
+// sparse_vector_type operator()(const sparse_vector_type& aa, const sparse_vector_type& bb);
+
+// value_type norm(const sparse_matrix_type& aa, Element_insertion_type it = Add_t);
+
+// value_type norm(const sparse_vector_type& aa, Element_insertion_type it = Add_t);
+
+// // NOTE: the destructor will return an error if it is called after MPI_Finalize is
+// // called because it uses collect communication to free-up allocated memory.
+// ~SolverPETSc() {
- static bool exit = false;
- if (!exit) {
- // add finalize PETSc function at exit
- atexit(finalize);
- exit = true;
- }
+// static bool exit = false;
+// if (!exit) {
+// // add finalize PETSc function at exit
+// atexit(finalize);
+// exit = true;
+// }
- if (allocated_) {
- PetscErrorCode ierr = MatDestroy(&A_);CHKERRCONTINUE(ierr);
- ierr = VecDestroy(&x_);CHKERRCONTINUE(ierr);
- ierr = KSPDestroy(&ksp_);CHKERRCONTINUE(ierr);
- }
- }
-
- /* from the PETSc library, these are the options that can be passed
- to the command line
+// if (allocated_) {
+// PetscErrorCode ierr = MatDestroy(&A_);CHKERRCONTINUE(ierr);
+// ierr = VecDestroy(&x_);CHKERRCONTINUE(ierr);
+// ierr = KSPDestroy(&ksp_);CHKERRCONTINUE(ierr);
+// }
+// }
+
+// /* from the PETSc library, these are the options that can be passed
+// to the command line
- Options Database Keys
+// Options Database Keys
- -options_table - Calls PetscOptionsView()
- -options_left - Prints unused options that remain in the database
- -objects_left - Prints list of all objects that have not been freed
- -mpidump - Calls PetscMPIDump()
- -malloc_dump - Calls PetscMallocDump()
- -malloc_info - Prints total memory usage
- -malloc_log - Prints summary of memory usage
+// -options_table - Calls PetscOptionsView()
+// -options_left - Prints unused options that remain in the database
+// -objects_left - Prints list of all objects that have not been freed
+// -mpidump - Calls PetscMPIDump()
+// -malloc_dump - Calls PetscMallocDump()
+// -malloc_info - Prints total memory usage
+// -malloc_log - Prints summary of memory usage
- Options Database Keys for Profiling
+// Options Database Keys for Profiling
- -log_summary [filename] - Prints summary of flop and timing information to screen.
- If the filename is specified the summary is written to the file. See PetscLogView().
- -log_summary_python [filename] - Prints data on of flop and timing usage to a file or screen.
- -log_all [filename] - Logs extensive profiling information See PetscLogDump().
- -log [filename] - Logs basic profiline information See PetscLogDump().
- -log_sync - Log the synchronization in scatters, inner products and norms
- -log_mpe [filename] - Creates a logfile viewable by the utility Upshot/Nupshot (in MPICH distribution)
- */
- static void finalize() {
+// -log_summary [filename] - Prints summary of flop and timing information to screen.
+// If the filename is specified the summary is written to the file. See PetscLogView().
+// -log_summary_python [filename] - Prints data on of flop and timing usage to a file or screen.
+// -log_all [filename] - Logs extensive profiling information See PetscLogDump().
+// -log [filename] - Logs basic profiline information See PetscLogDump().
+// -log_sync - Log the synchronization in scatters, inner products and norms
+// -log_mpe [filename] - Creates a logfile viewable by the utility Upshot/Nupshot (in MPICH distribution)
+// */
+// static void finalize() {
- static bool finalized = false;
- if (!finalized) {
+// static bool finalized = false;
+// if (!finalized) {
- cout<<"*** INFO *** PETSc is finalizing..."<<endl;
+// cout<<"*** INFO *** PETSc is finalizing..."<<endl;
- // finalize PETSc
- PetscErrorCode ierr = PetscFinalize();CHKERRCONTINUE(ierr);
- finalized = true;
+// // finalize PETSc
+// PetscErrorCode ierr = PetscFinalize();CHKERRCONTINUE(ierr);
+// finalized = true;
- cout<<"*** INFO *** Process "<<Parallel_base::rank_<<" is finalizing..."<<endl;
+// cout<<"*** INFO *** Process "<<Parallel_base::rank_<<" is finalizing..."<<endl;
- // finalize MPI
- MPI_Finalize();
- }
- }
-};
+// // finalize MPI
+// MPI_Finalize();
+// }
+// }
+// };
-class SolverPETSc : public Solver, public CommunicatorEventHandler {
- /* ------------------------------------------------------------------------ */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------ */
-public:
+// class SolverPETSc : public Solver, public CommunicatorEventHandler {
+// /* ------------------------------------------------------------------------ */
+// /* Constructors/Destructors */
+// /* ------------------------------------------------------------------------ */
+// public:
- SolverPETSc(SparseMatrix & sparse_matrix,
- const ID & id = "solver_petsc",
- const MemoryID & memory_id = 0);
+// SolverPETSc(SparseMatrix & sparse_matrix,
+// const ID & id = "solver_petsc",
+// const MemoryID & memory_id = 0);
- virtual SolverPETSc();
+// virtual SolverPETSc();
- /* ------------------------------------------------------------------------ */
- /* Methods */
- /* ------------------------------------------------------------------------ */
-public:
+// /* ------------------------------------------------------------------------ */
+// /* Methods */
+// /* ------------------------------------------------------------------------ */
+// public:
- /// build the profile and do the analysis part
- void initialize(SolverOptions & options = _solver_no_options);
+// /// build the profile and do the analysis part
+// void initialize(SolverOptions & options = _solver_no_options);
- void initializeSlave(SolverOptions & options = _solver_no_options);
+// void initializeSlave(SolverOptions & options = _solver_no_options);
- /// factorize and solve the system
- void solve(Array<Real> & solution);
- void solve();
+// /// factorize and solve the system
+// void solve(Array<Real> & solution);
+// void solve();
- void solveSlave();
+// void solveSlave();
- virtual void setRHS(Array<Real> & rhs);
+// virtual void setRHS(Array<Real> & rhs);
- /// function to print the contain of the class
- // virtual void printself(std::ostream & stream, int indent = 0) const;
+// /// function to print the contain of the class
+// // virtual void printself(std::ostream & stream, int indent = 0) const;
- virtual void onCommunicatorFinalize(const StaticCommunicator & communicator);
+// virtual void onCommunicatorFinalize(const StaticCommunicator & communicator);
-private:
+// private:
- void destroyMumpsData();
+// void destroyMumpsData();
- inline Int & icntl(UInt i) {
- return mumps_data.icntl[i - 1];
- }
+// inline Int & icntl(UInt i) {
+// return mumps_data.icntl[i - 1];
+// }
- inline Int & info(UInt i) {
- return mumps_data.info[i - 1];
- }
+// inline Int & info(UInt i) {
+// return mumps_data.info[i - 1];
+// }
- void initMumpsData(SolverMumpsOptions::ParallelMethod parallel_method);
+// void initMumpsData(SolverMumpsOptions::ParallelMethod parallel_method);
- /* ------------------------------------------------------------------------ */
- /* Accessors */
- /* ------------------------------------------------------------------------ */
-public:
+// /* ------------------------------------------------------------------------ */
+// /* Accessors */
+// /* ------------------------------------------------------------------------ */
+// public:
- /* ------------------------------------------------------------------------ */
- /* Class Members */
- /* ------------------------------------------------------------------------ */
-private:
+// /* ------------------------------------------------------------------------ */
+// /* Class Members */
+// /* ------------------------------------------------------------------------ */
+// private:
- /// mumps data
- DMUMPS_STRUC_C mumps_data;
+// /// mumps data
+// DMUMPS_STRUC_C mumps_data;
- /// specify if the mumps_data are initialized or not
- bool is_mumps_data_initialized;
+// /// specify if the mumps_data are initialized or not
+// bool is_mumps_data_initialized;
- UInt prank;
+// UInt prank;
- /* ------------------------------------------------------------------------ */
- /* Local types */
- /* ------------------------------------------------------------------------ */
+// /* ------------------------------------------------------------------------ */
+// /* Local types */
+// /* ------------------------------------------------------------------------ */
-private:
- SolverMumpsOptions::ParallelMethod parallel_method;
-
- bool rhs_is_local;
-
- enum SolverMumpsJob {
- _smj_initialize = -1,
- _smj_analyze = 1,
- _smj_factorize = 2,
- _smj_solve = 3,
- _smj_analyze_factorize = 4,
- _smj_factorize_solve = 5,
- _smj_complete = 6, // analyze, factorize, solve
- _smj_destroy = -2
- };
-};
+// private:
+// SolverMumpsOptions::ParallelMethod parallel_method;
+
+// bool rhs_is_local;
+
+// enum SolverMumpsJob {
+// _smj_initialize = -1,
+// _smj_analyze = 1,
+// _smj_factorize = 2,
+// _smj_solve = 3,
+// _smj_analyze_factorize = 4,
+// _smj_factorize_solve = 5,
+// _smj_complete = 6, // analyze, factorize, solve
+// _smj_destroy = -2
+// };
+// };
-/* -------------------------------------------------------------------------- */
-/* inline functions */
-/* -------------------------------------------------------------------------- */
+// /* -------------------------------------------------------------------------- */
+// /* inline functions */
+// /* -------------------------------------------------------------------------- */
-//#include "solver_mumps_inline_impl.cc"
+// //#include "solver_mumps_inline_impl.cc"
-/// standard output stream operator
-// inline std::ostream & operator <<(std::ostream & stream, const SolverMumps & _this)
-// {
-// _this.printself(stream);
-// return stream;
-// }
+// /// standard output stream operator
+// // inline std::ostream & operator <<(std::ostream & stream, const SolverMumps & _this)
+// // {
+// // _this.printself(stream);
+// // return stream;
+// // }
-__END_AKANTU__
+// __END_AKANTU__
-#endif /* __AKANTU_SOLVER_PETSC_HH__ */
+// #endif /* __AKANTU_SOLVER_PETSC_HH__ */
diff --git a/src/solver/sparse_matrix.hh b/src/solver/sparse_matrix.hh
index 6fed15019..21e06031a 100644
--- a/src/solver/sparse_matrix.hh
+++ b/src/solver/sparse_matrix.hh
@@ -1,253 +1,252 @@
/**
* @file sparse_matrix.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Mon Dec 13 2010
* @date last modification: Mon Sep 15 2014
*
* @brief sparse matrix storage class (distributed assembled matrix)
* This is a COO format (Coordinate List)
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_SPARSE_MATRIX_HH__
#define __AKANTU_SPARSE_MATRIX_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
#ifndef __INTEL_COMPILER
namespace std {
namespace tr1 {
template<typename a, typename b>
struct hash< std::pair<a, b> > {
private:
const hash<a> ah;
const hash<b> bh;
public:
hash() : ah(), bh() {}
size_t operator()(const std::pair<a, b> &p) const {
size_t seed = ah(p.first);
return bh(p.second) + 0x9e3779b9 + (seed<<6) + (seed>>2);
}
};
}
} // namespaces
#endif
__BEGIN_AKANTU__
class DOFSynchronizer;
class SparseMatrix : private Memory {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
SparseMatrix(UInt size,
const SparseMatrixType & sparse_matrix_type,
const ID & id = "sparse_matrix",
const MemoryID & memory_id = 0);
SparseMatrix(const SparseMatrix & matrix,
const ID & id = "sparse_matrix",
const MemoryID & memory_id = 0);
virtual ~SparseMatrix();
typedef std::pair<UInt, UInt> KeyCOO;
typedef unordered_map<KeyCOO, UInt>::type coordinate_list_map;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// remove the existing profile
inline void clearProfile();
/// add a non-zero element
- UInt addToProfile(UInt i, UInt j);
+ virtual UInt addToProfile(UInt i, UInt j);
/// set the matrix to 0
inline void clear();
/// assemble a local matrix in the sparse one
inline void addToMatrix(UInt i, UInt j, Real value);
/// set the size of the matrix
void resize(UInt size)
{ this->size = size; }
- /// fill the profil of the matrix
void buildProfile(const Mesh & mesh, const DOFSynchronizer & dof_synchronizer, UInt nb_degree_of_freedom);
/// modify the matrix to "remove" the blocked dof
- void applyBoundary(const Array<bool> & boundary, Real block_val = 1.);
+ virtual void applyBoundary(const Array<bool> & boundary, Real block_val = 1.);
// /// modify the matrix to "remove" the blocked dof
// void applyBoundaryNormal(Array<bool> & boundary_normal, Array<Real> & EulerAngles, Array<Real> & rhs, const Array<Real> & matrix, Array<Real> & rhs_rotated);
/// modify the matrix to "remove" the blocked dof
- void removeBoundary(const Array<bool> & boundary);
+ virtual void removeBoundary(const Array<bool> & boundary);
/// restore the profile that was before removing the boundaries
- void restoreProfile();
+ virtual void restoreProfile();
/// save the profil in a file using the MatrixMarket file format
- void saveProfile(const std::string & filename) const;
+ virtual void saveProfile(const std::string & filename) const;
/// save the matrix in a file using the MatrixMarket file format
- void saveMatrix(const std::string & filename) const;
+ virtual void saveMatrix(const std::string & filename) const;
/// copy assuming the profile are the same
- void copyContent(const SparseMatrix & matrix);
+ virtual void copyContent(const SparseMatrix & matrix);
/// copy profile
// void copyProfile(const SparseMatrix & matrix);
/// add matrix assuming the profile are the same
- void add(const SparseMatrix & matrix, Real alpha);
+ virtual void add(const SparseMatrix & matrix, Real alpha);
/// diagonal lumping
- void lump(Array<Real> & lumped);
+ virtual void lump(Array<Real> & lumped);
/// function to print the contain of the class
//virtual void printself(std::ostream & stream, int indent = 0) const;
-private:
+protected:
inline KeyCOO key(UInt i, UInt j) const {
if(sparse_matrix_type == _symmetric && (i > j))
return std::make_pair(j, i);
return std::make_pair(i, j);
}
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// return the values at potition i, j
inline Real operator()(UInt i, UInt j) const;
inline Real & operator()(UInt i, UInt j);
AKANTU_GET_MACRO(IRN, irn, const Array<Int> &);
AKANTU_GET_MACRO(JCN, jcn, const Array<Int> &);
AKANTU_GET_MACRO(A, a, const Array<Real> &);
AKANTU_GET_MACRO(NbNonZero, nb_non_zero, UInt);
AKANTU_GET_MACRO(Size, size, UInt);
AKANTU_GET_MACRO(SparseMatrixType, sparse_matrix_type, const SparseMatrixType &);
const DOFSynchronizer & getDOFSynchronizer() const {
AKANTU_DEBUG_ASSERT(dof_synchronizer != NULL,
"DOFSynchronizer not initialized in the SparseMatrix!");
return *dof_synchronizer;
}
private:
AKANTU_GET_MACRO(DOFSynchronizerPointer, dof_synchronizer, DOFSynchronizer *);
friend Array<Real> & operator*=(Array<Real> & vect, const SparseMatrix & mat);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
-private:
+protected:
/// id of the SparseMatrix
ID id;
/// sparce matrix type
SparseMatrixType sparse_matrix_type;
/// Mesh corresponding to the profile
// const Mesh * mesh;
/// Size of the matrix
UInt size;
/// number of processors
UInt nb_proc;
/// number of non zero element
UInt nb_non_zero;
/// row indexes
Array<Int> irn;
/// column indexes
Array<Int> jcn;
/// values : A[k] = Matrix[irn[k]][jcn[k]]
Array<Real> a;
/// saved row indexes
Array<Int> * irn_save;
/// saved column indexes
Array<Int> * jcn_save;
/// saved size
UInt size_save;
/// information to know where to assemble an element in a global sparse matrix
// ElementTypeMapArray<UInt> element_to_sparse_profile;
/* map for (i,j) -> k correspondence \warning std::map are slow
* \todo improve with hash_map (non standard in stl) or unordered_map (boost or C++0x)
*/
coordinate_list_map irn_jcn_k;
DOFSynchronizer * dof_synchronizer;
// std::map<std::pair<UInt, UInt>, UInt> * irn_jcn_to_k;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "sparse_matrix_inline_impl.cc"
#endif
// /// standard output stream operator
// inline std::ostream & operator <<(std::ostream & stream, const SparseMatrix & _this)
// {
// _this.printself(stream);
// return stream;
// }
Array<Real> & operator*=(Array<Real> & vect, const SparseMatrix & mat);
__END_AKANTU__
#endif /* __AKANTU_SPARSE_MATRIX_HH__ */
diff --git a/src/solver/static_solver.cc b/src/solver/static_solver.cc
new file mode 100644
index 000000000..be55b89b2
--- /dev/null
+++ b/src/solver/static_solver.cc
@@ -0,0 +1,102 @@
+/**
+ * @file static_solver.cc
+ * @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
+ * @date Wed Jul 30 15:35:01 2014
+ *
+ * @brief implementation of the static solver
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "static_solver.hh"
+/* -------------------------------------------------------------------------- */
+#ifdef AKANTU_USE_PETSC
+#include <petscsys.h>
+#endif
+
+/* -------------------------------------------------------------------------- */
+
+__BEGIN_AKANTU__
+
+/* -------------------------------------------------------------------------- */
+StaticSolver::~StaticSolver() {
+ --this->nb_references;
+ if(this->nb_references == 0)
+ delete this->static_solver;
+}
+
+/* -------------------------------------------------------------------------- */
+StaticSolver & StaticSolver::getStaticSolver() {
+ if(nb_references == 0)
+ static_solver = new StaticSolver();
+
+ ++nb_references;
+
+ return *static_solver;
+}
+
+#ifdef AKANTU_USE_PETSC
+#if PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 5
+static PetscErrorCode PETScErrorHandler(MPI_Comm,
+ int line, const char * dir, const char *file,
+ PetscErrorCode number,
+ PetscErrorType type,
+ const char *message,
+ void *) {
+ AKANTU_DEBUG_ERROR("An error occured in PETSc in file \"" << file << ":" << line << "\" - PetscErrorCode "<< number << " - \""<< message << "\"");
+}
+#else
+static PetscErrorCode PETScErrorHandler(MPI_Comm,
+ int line, const char * func, const char * dir, const char *file,
+ PetscErrorCode number,
+ PetscErrorType type,
+ const char *message,
+ void *) {
+ AKANTU_DEBUG_ERROR("An error occured in PETSc in file \"" << file << ":" << line << "\" - PetscErrorCode "<< number << " - \""<< message << "\"");
+}
+#endif
+#endif
+
+/* -------------------------------------------------------------------------- */
+void StaticSolver::initialize(int & argc, char ** & argv) {
+ AKANTU_DEBUG_ASSERT(this->is_initialized != true, "The static solver has already been initialized");
+#ifdef AKANTU_USE_PETSC
+ PetscErrorCode petsc_error = PetscInitialize(&argc, &argv, NULL, NULL);
+ if(petsc_error != 0) {
+ AKANTU_DEBUG_ERROR("An error occured while initializing Petsc (PetscErrorCode "<< petsc_error << ")");
+ }
+ PetscPushErrorHandler(PETScErrorHandler, NULL);
+#endif
+
+ this->is_initialized = true;
+ }
+
+/* -------------------------------------------------------------------------- */
+void StaticSolver::finalize() {
+ AKANTU_DEBUG_ASSERT(this->is_initialized == true, "The static solver has not been initialized");
+#ifdef AKANTU_USE_PETSC
+ PetscFinalize();
+#endif
+
+ this->is_initialized = false;
+}
+
+__END_AKANTU__
diff --git a/src/solver/static_solver.hh b/src/solver/static_solver.hh
new file mode 100644
index 000000000..add789674
--- /dev/null
+++ b/src/solver/static_solver.hh
@@ -0,0 +1,73 @@
+/**
+ * @file static_solver.hh
+ * @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
+ * @date Wed Jul 30 14:44:12 2014
+ *
+ * @brief Class handeling the initialization of external solvers
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "aka_common.hh"
+
+
+#ifndef __AKANTU_STATIC_SOLVER_HH__
+#define __AKANTU_STATIC_SOLVER_HH__
+
+__BEGIN_AKANTU__
+
+class StaticSolver {
+ /* ------------------------------------------------------------------------ */
+ /* Constructors */
+ /* ------------------------------------------------------------------------ */
+private:
+ StaticSolver() : is_initialized(false) {};
+
+public:
+ ~StaticSolver();
+
+ /* ------------------------------------------------------------------------ */
+ /// get an instance to the static solver
+ static StaticSolver & getStaticSolver();
+
+ /* ------------------------------------------------------------------------ */
+ /* Methods */
+ /* ------------------------------------------------------------------------ */
+public:
+ /// initialize what is needed for the compiled solver interfaces
+ void initialize(int & argc, char ** & argv);
+
+ /// finalize what is needed for the compiled solver interfaces
+ void finalize();
+
+ /* ------------------------------------------------------------------------ */
+ /* Members */
+ /* ------------------------------------------------------------------------ */
+private:
+ bool is_initialized;
+ static UInt nb_references;
+ static StaticSolver * static_solver;
+};
+
+__END_AKANTU__
+
+
+#endif /* __AKANTU_STATIC_SOLVER_HH__ */
diff --git a/test/test_solver/CMakeLists.txt b/test/test_solver/CMakeLists.txt
index 9f4b59df5..42f17633c 100644
--- a/test/test_solver/CMakeLists.txt
+++ b/test/test_solver/CMakeLists.txt
@@ -1,49 +1,54 @@
#===============================================================================
# @file CMakeLists.txt
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date creation: Mon Dec 13 2010
# @date last modification: Tue Nov 06 2012
#
# @brief configuration for solver tests
#
# @section LICENSE
#
# Copyright (©) 2014 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 <http://www.gnu.org/licenses/>.
#
# @section DESCRIPTION
#
#===============================================================================
add_mesh(test_solver_mesh triangle.geo 2 1)
register_test(test_sparse_matrix_profile
SOURCES test_sparse_matrix_profile.cc
DEPENDENCIES test_solver_mesh
)
register_test(test_sparse_matrix_assemble
SOURCES test_sparse_matrix_assemble.cc
DEPENDENCIES test_solver_mesh
)
register_test(test_sparse_matrix_product
SOURCES test_sparse_matrix_product.cc
FILES_TO_COPY bar.msh
)
+
+register_test(test_petsc_matrix_profile
+ SOURCES test_petsc_matrix_profile.cc
+ DEPENDENCIES test_solver_mesh
+ )
\ No newline at end of file
diff --git a/test/test_solver/test_petsc_matrix_profile.cc b/test/test_solver/test_petsc_matrix_profile.cc
new file mode 100644
index 000000000..22c0e90d6
--- /dev/null
+++ b/test/test_solver/test_petsc_matrix_profile.cc
@@ -0,0 +1,75 @@
+/**
+ * @file test_petsc_matrix_profile.cc
+ * @author Aurelia Cuba Ramos <aurelia.cubaramos@epfl.ch>
+ * @date Wed Jul 30 12:34:08 2014
+ *
+ * @brief test the profile generation of the PETScMatrix class
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2010-2011 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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include <cstdlib>
+/* -------------------------------------------------------------------------- */
+#include "static_communicator.hh"
+#include "aka_common.hh"
+#include "mesh.hh"
+#include "mesh_io.hh"
+
+#include "petsc_matrix.hh"
+#include "dof_synchronizer.hh"
+
+#include "mesh_partition_scotch.hh"
+
+using namespace akantu;
+int main(int argc, char *argv[]) {
+
+ initialize(argc, argv);
+ UInt spatial_dimension = 2;
+
+ StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator();
+ Int psize = comm.getNbProc();
+ Int prank = comm.whoAmI();
+
+ /// read the mesh and partition it
+ Mesh mesh(spatial_dimension);
+ akantu::MeshPartition * partition = NULL;
+
+ if(prank == 0) {
+ /// creation mesh
+ mesh.read("triangle.msh");
+ partition = new MeshPartitionScotch(mesh, spatial_dimension);
+ partition->partitionate(psize);
+ }
+
+
+ DOFSynchronizer dof_synchronizer(mesh, spatial_dimension);
+ UInt nb_global_nodes = mesh.getNbGlobalNodes();
+
+ PETScMatrix petsc_matrix(nb_global_nodes * spatial_dimension, _symmetric);
+
+ dof_synchronizer.initGlobalDOFEquationNumbers();
+
+ petsc_matrix.resize(dof_synchronizer);
+ petsc_matrix.buildProfile(mesh, dof_synchronizer, spatial_dimension);
+
+ petsc_matrix.saveMatrix("profile.mtx");
+
+}