diff --git a/src/common/aka_array.hh b/src/common/aka_array.hh index a77eaf50e..5d5fce77e 100644 --- a/src/common/aka_array.hh +++ b/src/common/aka_array.hh @@ -1,373 +1,375 @@ /** * @file aka_array.hh * * @author Till Junge * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Jun 24 2014 * * @brief Array container for Akantu * This container differs from the std::vector from the fact it as 2 dimensions * a main dimension and the size stored per entries * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_VECTOR_HH__ #define __AKANTU_VECTOR_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #include //#include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /// class that afford to store vectors in static memory class ArrayBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ArrayBase(const ID & id = ""); virtual ~ArrayBase(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get the amount of space allocated in bytes inline UInt getMemorySize() const; /// set the size to zero without freeing the allocated space inline void empty(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get the real size allocated in memory AKANTU_GET_MACRO(AllocatedSize, allocated_size, UInt); /// Get the Size of the Array AKANTU_GET_MACRO(Size, size, UInt); /// Get the number of components AKANTU_GET_MACRO(NbComponent, nb_component, UInt); /// Get the name of th array AKANTU_GET_MACRO(ID, id, const ID &); + /// Set the name of th array + AKANTU_SET_MACRO(ID, id, const ID &); // AKANTU_GET_MACRO(Tag, tag, const std::string &); // AKANTU_SET_MACRO(Tag, tag, const std::string &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the vector ID id; /// the size allocated UInt allocated_size; /// the size used UInt size; /// number of components UInt nb_component; /// size of the stored type UInt size_of_type; // /// User defined tag // std::string tag; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class Array : public ArrayBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef T value_type; typedef value_type & reference; typedef value_type * pointer_type; typedef const value_type & const_reference; /// Allocation of a new vector inline Array(UInt size = 0, UInt nb_component = 1, const ID & id = ""); /// Allocation of a new vector with a default value Array(UInt size, UInt nb_component, const value_type def_values[], const ID & id = ""); /// Allocation of a new vector with a default value Array(UInt size, UInt nb_component, const_reference value, const ID & id = ""); /// Copy constructor (deep copy if deep=true) Array(const Array& vect, bool deep = true, const ID & id = ""); #ifndef SWIG /// Copy constructor (deep copy) Array(const std::vector & vect); #endif virtual inline ~Array(); Array & operator=(const Array & a) { /// this is to let STL allocate and copy arrays in the case of std::vector::resize AKANTU_DEBUG_ASSERT(this->size == 0,"Cannot copy akantu::Array"); return const_cast(a); } /* ------------------------------------------------------------------------ */ /* Iterator */ /* ------------------------------------------------------------------------ */ /// \todo protected: does not compile with intel check why public: template ::value > class iterator_internal; public: /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ template class const_iterator; template class iterator; /* ------------------------------------------------------------------------ */ /// iterator for Array of nb_component = 1 typedef iterator< T > scalar_iterator; /// const_iterator for Array of nb_component = 1 typedef const_iterator< T > const_scalar_iterator; /// iterator rerturning Vectors of size n on entries of Array with nb_component = n typedef iterator< Vector > vector_iterator; /// const_iterator rerturning Vectors of n size on entries of Array with nb_component = n typedef const_iterator< Vector > const_vector_iterator; /// iterator rerturning Matrices of size (m, n) on entries of Array with nb_component = m*n typedef iterator< Matrix > matrix_iterator; /// const iterator rerturning Matrices of size (m, n) on entries of Array with nb_component = m*n typedef const_iterator< Matrix > const_matrix_iterator; /* ------------------------------------------------------------------------ */ /// Get an iterator that behaves like a pointer T * to the first entry inline iterator begin(); /// Get an iterator that behaves like a pointer T * to the end of the Array inline iterator end(); /// Get a const_iterator to the beginging of an Array of scalar inline const_iterator begin() const; /// Get a const_iterator to the end of an Array of scalar inline const_iterator end() const; /* ------------------------------------------------------------------------ */ /// Get a vector_iterator on the begining of the Array inline vector_iterator begin(UInt n); /// Get a vector_iterator on the end of the Array inline vector_iterator end(UInt n); /// Get a vector_iterator on the begining of the Array inline const_vector_iterator begin(UInt n) const; /// Get a vector_iterator on the end of the Array inline const_vector_iterator end(UInt n) const; /// Get a vector_iterator on the begining of the Array considered of shape (new_size, n) inline vector_iterator begin_reinterpret(UInt n, UInt new_size); /// Get a vector_iterator on the end of the Array considered of shape (new_size, n) inline vector_iterator end_reinterpret(UInt n, UInt new_size); /// Get a const_vector_iterator on the begining of the Array considered of shape (new_size, n) inline const_vector_iterator begin_reinterpret(UInt n, UInt new_size) const; /// Get a const_vector_iterator on the end of the Array considered of shape (new_size, n) inline const_vector_iterator end_reinterpret(UInt n, UInt new_size) const; /* ------------------------------------------------------------------------ */ /// Get a matrix_iterator on the begining of the Array (Matrices of size (m, n)) inline matrix_iterator begin(UInt m, UInt n); /// Get a matrix_iterator on the end of the Array (Matrices of size (m, n)) inline matrix_iterator end(UInt m, UInt n); /// Get a const_matrix_iterator on the begining of the Array (Matrices of size (m, n)) inline const_matrix_iterator begin(UInt m, UInt n) const; /// Get a const_matrix_iterator on the end of the Array (Matrices of size (m, n)) inline const_matrix_iterator end(UInt m, UInt n) const; /// Get a matrix_iterator on the begining of the Array considered of shape (new_size, m*n) inline matrix_iterator begin_reinterpret(UInt m, UInt n, UInt size); /// Get a matrix_iterator on the end of the Array considered of shape (new_size, m*n) inline matrix_iterator end_reinterpret(UInt m, UInt n, UInt size); /// Get a const_matrix_iterator on the begining of the Array considered of shape (new_size, m*n) inline const_matrix_iterator begin_reinterpret(UInt m, UInt n, UInt size) const; /// Get a const_matrix_iterator on the end of the Array considered of shape (new_size, m*n) inline const_matrix_iterator end_reinterpret(UInt m, UInt n, UInt size) const; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// append a tuple of size nb_component containing value inline void push_back(const_reference value); /// append a vector inline void push_back(const value_type new_elem[]); /// append a Vector or a Matrix template class C> inline void push_back(const C & new_elem); /// append the value of the iterator template inline void push_back(const iterator & it); /// erase the value at position i inline void erase(UInt i); /// ask Nico, clarify template inline iterator erase(const iterator & it); /// change the size of the Array void resize(UInt size); /// change the number of components by interlacing data void extendComponentsInterlaced(UInt multiplicator, UInt stride); /// search elem in the vector, return the position of the first occurrence or -1 if not found Int find(const_reference elem) const;\ /// @see Array::find(const_reference elem) const Int find(T elem[]) const; /// set all entries of the array to 0 inline void clear() { std::fill_n(values, size*nb_component, T()); } /// set all entries of the array to the value t /// @param t value to fill the array with inline void set(T t) { std::fill_n(values, size*nb_component, t); } /// set all tuples of the array to a given vector or matrix /// @param vm Matrix or Vector to fill the array with template class C> inline void set(const C & vm); /// copy another Array in the current Array, the no_sanity_check allows you to /// force the copy in cases where you know what you do with two non matching /// Arrays in terms of n void copy(const Array & other, bool no_sanity_check = false); /// give the address of the memory allocated for this vector T * storage() const { return values; }; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// perform the allocation for the constructors void allocate(UInt size, UInt nb_component = 1); /// resize without initializing the memory void resizeUnitialized(UInt new_size); /* ------------------------------------------------------------------------ */ /* Operators */ /* ------------------------------------------------------------------------ */ public: /// substraction entry-wise Array & operator-=(const Array & other); /// addition entry-wise Array & operator+=(const Array & other); /// multiply evry entry by alpha Array & operator*=(const T & alpha); /// check if the array are identical entry-wise bool operator==(const Array & other) const; /// @see Array::operator==(const Array & other) const bool operator!=(const Array & other) const; /// return a reference to the j-th entry of the i-th tuple inline reference operator()(UInt i, UInt j = 0); /// return a const reference to the j-th entry of the i-th tuple inline const_reference operator()(UInt i, UInt j = 0) const; /// return a reference to the ith component of the 1D array inline reference operator[](UInt i); /// return a const reference to the ith component of the 1D array inline const_reference operator[](UInt i) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the number of tuples contained in the array UInt getSize() const{ return this->size; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// array of values T * values; // /!\ very dangerous }; __END_AKANTU__ #include "aka_types.hh" __BEGIN_AKANTU__ #include "aka_array_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Inline Functions Array */ /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Array & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ /* Inline Functions ArrayBase */ /* -------------------------------------------------------------------------- */ inline std::ostream & operator<<(std::ostream & stream, const ArrayBase & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_VECTOR_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc index 5b4c0b21a..f5673a24c 100644 --- a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc @@ -1,823 +1,832 @@ /** * @file material_non_local_inline_impl.cc * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Mon Jun 23 2014 * * @brief Non-local inline implementation * * @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 . * */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ /* -------------------------------------------------------------------------- */ #include "aka_types.hh" #include "grid_synchronizer.hh" #include "synchronizer_registry.hh" #include "integrator.hh" +#include "dumper_paraview.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #if defined(AKANTU_DEBUG_TOOLS) # include "aka_debug_tools.hh" #endif __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::MaterialNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), weight_func(NULL), spatial_grid(NULL), compute_stress_calls(0), is_creating_grid(false), grid_synchronizer(NULL) { AKANTU_DEBUG_IN(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; pair_weight[ghost_type] = NULL; } this->is_non_local = true; this->weight_func = new WeightFunction(*this); this->registerSubSection(_st_non_local, "weight_function", *weight_func); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::~MaterialNonLocal() { AKANTU_DEBUG_IN(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; delete pair_weight[ghost_type]; } delete spatial_grid; delete weight_func; delete grid_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::initMaterial() { AKANTU_DEBUG_IN(); // Material::initMaterial(); Mesh & mesh = this->model->getFEEngine().getMesh(); InternalField quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", *this); quadrature_points_coordinates.initialize(spatial_dimension); ElementTypeMap nb_ghost_protected; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); for(; it != last_type; ++it) nb_ghost_protected(mesh.getNbElement(*it, _ghost), *it, _ghost); AKANTU_DEBUG_INFO("Creating cell list"); createCellList(quadrature_points_coordinates); AKANTU_DEBUG_INFO("Creating pairs"); updatePairList(quadrature_points_coordinates); #if not defined(AKANTU_NDEBUG) if(AKANTU_DEBUG_TEST(dblDump)) neighbourhoodStatistics("material_non_local.stats"); #endif AKANTU_DEBUG_INFO("Cleaning extra ghosts"); - cleanupExtraGhostElement(nb_ghost_protected); + /// cleanupExtraGhostElement(nb_ghost_protected); AKANTU_DEBUG_INFO("Computing weights"); weight_func->setRadius(weight_func->getRadius()); weight_func->init(); computeWeights(quadrature_points_coordinates); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::cleanupExtraGhostElement(const ElementTypeMap & nb_ghost_protected) { AKANTU_DEBUG_IN(); // Create list of element to keep std::set relevant_ghost_element; PairList::const_iterator first_pair = pair_list[_ghost].begin(); PairList::const_iterator last_pair = pair_list[_ghost].end(); for(;first_pair != last_pair; ++first_pair) { const QuadraturePoint & q2 = first_pair->second; relevant_ghost_element.insert(q2); } // Create list of element to remove and new numbering for element to keep Mesh & mesh = this->model->getFEEngine().getMesh(); std::set ghost_to_erase; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); RemovedElementsEvent remove_elem(mesh); - Element element_global; // member element corresponds to global element number - element_global.ghost_type = _ghost; + Element element_local; // member element corresponds to global element number + element_local.ghost_type = _ghost; for(; it != last_type; ++it) { - element_global.type = *it; - UInt nb_ghost_elem = mesh.getNbElement(*it, _ghost); + element_local.type = *it; + UInt nb_ghost_elem_global = mesh.getNbElement(*it, _ghost); UInt nb_ghost_elem_protected = 0; try { nb_ghost_elem_protected = nb_ghost_protected(*it, _ghost); } catch (...) {} if(!remove_elem.getNewNumbering().exists(*it, _ghost)) - remove_elem.getNewNumbering().alloc(nb_ghost_elem, 1, *it, _ghost); - else remove_elem.getNewNumbering(*it, _ghost).resize(nb_ghost_elem); - + remove_elem.getNewNumbering().alloc(nb_ghost_elem_global, 1, *it, _ghost); + else remove_elem.getNewNumbering(*it, _ghost).resize(nb_ghost_elem_global); + Array & elem_filter = element_filter(*it, _ghost); + UInt nb_ghost_elem_local = elem_filter.getSize(); Array & new_numbering = remove_elem.getNewNumbering(*it, _ghost); - UInt ng = 0; - for (UInt g = 0; g < nb_ghost_elem; ++g) { - if (element_global.element >= nb_ghost_elem_protected) { - Element element_local = this->convertToLocalElement(element_global); - - if (std::find(relevant_ghost_element.begin(), - relevant_ghost_element.end(), - element_local) == relevant_ghost_element.end()) { + for (UInt g = 0; g < nb_ghost_elem_local; ++g) { + element_local.element = g; + Element element_global = this->convertToGlobalElement(element_local); + if (element_global.element >= nb_ghost_elem_protected && + relevant_ghost_element.find(element_local) == relevant_ghost_element.end()) { + // (std::find(relevant_ghost_element.begin(), + // relevant_ghost_element.end(), + // element_local) == relevant_ghost_element.end())) { + std::cout<< "element removed" << std::endl; remove_elem.getList().push_back(element_global); - new_numbering(g) = UInt(-1); + new_numbering(element_global.element) = UInt(-1); } - } else { - new_numbering(g) = ng; + } + + UInt ng = 0; + for (UInt g = 0; g < nb_ghost_elem_global; ++g) { + if (new_numbering(g) != UInt(-1)) { + new_numbering(g) = ng; ++ng; } } } mesh.sendEvent(remove_elem); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::createCellList(ElementTypeMapArray & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); const Real safety_factor = 1.2; // for the cell grid spacing Mesh & mesh = this->model->getFEEngine().getMesh(); mesh.computeBoundingBox(); const Vector & lower_bounds = mesh.getLocalLowerBounds(); const Vector & upper_bounds = mesh.getLocalUpperBounds(); Vector center = 0.5 * (upper_bounds + lower_bounds); Vector spacing(spatial_dimension, weight_func->getRadius() * safety_factor); spatial_grid = new SpatialGrid(spatial_dimension, spacing, center); this->computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); this->fillCellList(quadrature_points_coordinates, _not_ghost); - + DumperParaview dumper_ghost("ghosts"); + dumper_ghost.registerMesh(mesh, spatial_dimension, _ghost); + dumper_ghost.dump(); is_creating_grid = true; std::set tags; tags.insert(_gst_mnl_for_average); tags.insert(_gst_mnl_weight); tags.insert(_gst_material_id); SynchronizerRegistry & synch_registry = this->model->getSynchronizerRegistry(); std::stringstream sstr; sstr << getID() << ":grid_synchronizer"; grid_synchronizer = GridSynchronizer::createGridSynchronizer(mesh, *spatial_grid, sstr.str(), &synch_registry, tags); is_creating_grid = false; this->computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); fillCellList(quadrature_points_coordinates, _ghost); - + dumper_ghost.dump(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::fillCellList(const ElementTypeMapArray & quadrature_points_coordinates, const GhostType & ghost_type) { QuadraturePoint q; q.ghost_type = ghost_type; Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->element_filter.lastType (spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quad = this->model->getFEEngine().getNbQuadraturePoints(*it, ghost_type); const Array & quads = quadrature_points_coordinates(*it, ghost_type); q.type = *it; Array::const_vector_iterator quad = quads.begin(spatial_dimension); UInt * elem = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e) { q.element = *elem; for (UInt nq = 0; nq < nb_quad; ++nq) { q.num_point = nq; q.global_num = q.element * nb_quad + nq; //q.setPosition(*quad); spatial_grid->insert(q, *quad); ++quad; } ++elem; } } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updatePairList(const ElementTypeMapArray & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); GhostType ghost_type = _not_ghost; QuadraturePoint q1; q1.ghost_type = ghost_type; // generate the pair of neighbor depending of the cell_list Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { // Preparing datas const Array & quads = quadrature_points_coordinates(*it, ghost_type); Array::const_vector_iterator q1_coord_it = quads.begin(spatial_dimension); Array::const_vector_iterator last_quad = quads.end(spatial_dimension); q1.type = *it; q1.global_num = 0; // loop over quad points for(;q1_coord_it != last_quad; ++q1_coord_it, ++(q1.global_num)) { UInt nb_quad1 = this->model->getFEEngine().getNbQuadraturePoints(q1.type, q1.ghost_type); q1.element = q1.global_num / nb_quad1; q1.num_point = q1.global_num % nb_quad1; const Vector & q1_coord = *q1_coord_it; SpatialGrid::CellID cell_id = spatial_grid->getCellID(q1_coord); SpatialGrid::neighbor_cells_iterator first_neigh_cell = spatial_grid->beginNeighborCells(cell_id); SpatialGrid::neighbor_cells_iterator last_neigh_cell = spatial_grid->endNeighborCells(cell_id); // loop over neighbors cells of the one containing the current quadrature // point for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { SpatialGrid::Cell::iterator first_neigh_quad = spatial_grid->beginCell(*first_neigh_cell); SpatialGrid::Cell::iterator last_neigh_quad = spatial_grid->endCell(*first_neigh_cell); // loop over the quadrature point in the current cell of the cell list for (;first_neigh_quad != last_neigh_quad; ++first_neigh_quad){ QuadraturePoint q2 = this->convertToLocalPoint(*first_neigh_quad); Array::const_vector_iterator q2_coord_it = quadrature_points_coordinates(q2.type, q2.ghost_type).begin(spatial_dimension); const Vector & q2_coord = q2_coord_it[q2.global_num]; Real distance = q1_coord.distance(q2_coord); if(distance <= weight_func->getRadius() && (q2.ghost_type == _ghost || (q2.ghost_type == _not_ghost && q1.global_num <= q2.global_num))) { // storing only half lists pair_list[q2.ghost_type].push_back(std::make_pair(q1, q2)); } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::computeWeights(const ElementTypeMapArray & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); InternalField quadrature_points_volumes("quadrature_points_volumes", *this); quadrature_points_volumes.initialize(1); const FEEngine & fem = this->model->getFEEngine(); weight_func->updateInternals(quadrature_points_volumes); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; if(!(pair_weight[ghost_type2])) { std::string ghost_id = ""; if (ghost_type2 == _ghost) ghost_id = ":ghost"; std::stringstream sstr; sstr << getID() << ":pair_weight:" << ghost_id; pair_weight[ghost_type2] = new Array(0, 2, sstr.str()); } pair_weight[ghost_type2]->resize(pair_list[ghost_type2].size()); pair_weight[ghost_type2]->clear(); PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); // Compute the weights for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & lq1 = first_pair->first; const QuadraturePoint & lq2 = first_pair->second; QuadraturePoint gq1 = this->convertToGlobalPoint(lq1); QuadraturePoint gq2 = this->convertToGlobalPoint(lq2); // const Real q2_wJ = fem.getIntegratorInterface().getJacobians(gq2.type, gq2.ghost_type)(gq2.global_num); Array::const_vector_iterator quad_coords_1 = quadrature_points_coordinates(lq1.type, lq1.ghost_type).begin(spatial_dimension); Array::const_vector_iterator quad_coords_2 = quadrature_points_coordinates(lq2.type, lq2.ghost_type).begin(spatial_dimension); Array & quad_volumes_1 = quadrature_points_volumes(lq1.type, lq1.ghost_type); const Array & jacobians_2 = fem.getIntegratorInterface().getJacobians(gq2.type, gq2.ghost_type); const Real & q2_wJ = jacobians_2(gq2.global_num); // Real & q1_volume = quad_volumes(lq1.global_num); const Vector & q1_coord = quad_coords_1[lq1.global_num]; // quadrature_points_coordinates(lq1.type, lq1.ghost_type).begin(spatial_dimension)[lq1.global_num]; const Vector & q2_coord = quad_coords_2[lq2.global_num]; // quadrature_points_coordinates(lq1.type, lq1.ghost_type).begin(spatial_dimension)[lq2.global_num]; this->weight_func->selectType(lq1.type, lq1.ghost_type, lq2.type, lq2.ghost_type); // Weight function Real r = q1_coord.distance(q2_coord); Real w1 = this->weight_func->operator()(r, lq1, lq2); weight(0) = q2_wJ * w1; // q1_volume += weight(0); quad_volumes_1(lq1.global_num) += weight(0); if(lq2.ghost_type != _ghost && lq1.global_num != lq2.global_num) { const Array & jacobians_1 = fem.getIntegratorInterface().getJacobians(gq1.type, gq1.ghost_type); Array & quad_volumes_2 = quadrature_points_volumes(lq2.type, lq2.ghost_type); const Real & q1_wJ = jacobians_1(gq1.global_num); //Real & q2_volume = quad_volumes_2(lq2.global_num); Real w2 = this->weight_func->operator()(r, lq2, lq1); weight(1) = q1_wJ * w2; quad_volumes_2(lq2.global_num) += weight(1); //q2_volume += weight(1); } else weight(1) = 0.; } } //normalize the weights for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); // Compute the weights for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & lq1 = first_pair->first; const QuadraturePoint & lq2 = first_pair->second; Array & quad_volumes_1 = quadrature_points_volumes(lq1.type, lq1.ghost_type); Array & quad_volumes_2 = quadrature_points_volumes(lq2.type, lq2.ghost_type); Real q1_volume = quad_volumes_1(lq1.global_num); weight(0) *= 1. / q1_volume; if(ghost_type2 != _ghost) { Real q2_volume = quad_volumes_2(lq2.global_num); weight(1) *= 1. / q2_volume; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> template void MaterialNonLocal::weightedAvergageOnNeighbours(const InternalField & to_accumulate, InternalField & accumulated, UInt nb_degree_of_freedom, GhostType ghost_type2) const { AKANTU_DEBUG_IN(); if(ghost_type2 == _not_ghost) { accumulated.reset(); } PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); // Compute the weights for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & lq1 = first_pair->first; const QuadraturePoint & lq2 = first_pair->second; const Array & to_acc_1 = to_accumulate(lq1.type, lq1.ghost_type); Array & acc_1 = accumulated(lq1.type, lq1.ghost_type); const Array & to_acc_2 = to_accumulate(lq2.type, lq2.ghost_type); Array & acc_2 = accumulated(lq2.type, lq2.ghost_type); // const Vector & q2_to_acc = to_acc_2[lq2.global_num]; // Vector & q1_acc = acc_1[lq1.global_num]; //q1_acc += weight(0) * q2_to_acc; for(UInt d = 0; d < nb_degree_of_freedom; ++d) { acc_1(lq1.global_num, d) += weight(0) * to_acc_2(lq2.global_num, d); } if(ghost_type2 != _ghost) { for(UInt d = 0; d < nb_degree_of_freedom; ++d) { acc_2(lq2.global_num, d) += weight(1) * to_acc_1(lq1.global_num, d); } } // if(ghost_type2 != _ghost) { // const Vector & q1_to_acc = to_acc_1[lq1.global_num]; // Vector & q2_acc = acc_2[lq2.global_num]; // q2_acc += weight(1) * q1_to_acc; // } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); // Update the weights for the non local variable averaging if(ghost_type == _not_ghost && this->weight_func->getUpdateRate() && (this->compute_stress_calls % this->weight_func->getUpdateRate() == 0)) { ElementTypeMapArray quadrature_points_coordinates("quadrature_points_coordinates", getID()); Mesh & mesh = this->model->getFEEngine().getMesh(); mesh.initElementTypeMapArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); computeWeights(quadrature_points_coordinates); } if(ghost_type == _not_ghost) ++this->compute_stress_calls; computeAllStresses(ghost_type); computeNonLocalStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::computeAllNonLocalStresses(GhostType ghost_type) { // Update the weights for the non local variable averaging if(ghost_type == _not_ghost) { if(this->weight_func->getUpdateRate() && (this->compute_stress_calls % this->weight_func->getUpdateRate() == 0)) { this->model->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_weight); ElementTypeMapArray quadrature_points_coordinates("quadrature_points_coordinates", getID()); Mesh & mesh = this->model->getFEEngine().getMesh(); mesh.initElementTypeMapArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); this->model->getSynchronizerRegistry().waitEndSynchronize(_gst_mnl_weight); computeWeights(quadrature_points_coordinates); } typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; non_local_variable.non_local->resize(); this->weightedAvergageOnNeighbours(*non_local_variable.local, *non_local_variable.non_local, non_local_variable.nb_component, _not_ghost); } ++this->compute_stress_calls; } else { typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; this->weightedAvergageOnNeighbours(*non_local_variable.local, *non_local_variable.non_local, non_local_variable.nb_component, _ghost); } computeNonLocalStresses(_not_ghost); } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::savePairs(const std::string & filename) const { std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); for(;first_pair != last_pair; ++first_pair, ++weight_it) { Vector & weight = *weight_it; const QuadraturePoint & lq1 = first_pair->first; const QuadraturePoint & lq2 = first_pair->second; pout << lq1 << " " << lq2 << " " << weight << std::endl; } } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::neighbourhoodStatistics(const std::string & filename) const { // std::ofstream pout; // pout.open(filename.c_str()); // const Mesh & mesh = this->model->getFEEngine().getMesh(); // GhostType ghost_type1; // ghost_type1 = _not_ghost; // StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); // Int prank = comm.whoAmI(); // InternalField nb_neighbors("nb_neighbours", *const_cast(this)); // nb_neighbors.initialize(1); // for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { // GhostType ghost_type2 = (GhostType) gt; // UInt existing_pairs_num = gt - _not_ghost; // pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); // pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); // for (; first_pair_types != last_pair_types; ++first_pair_types) { // const Array & pairs = // pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); // if(prank == 0) { // pout << ghost_type2 << " "; // pout << "Types : " << first_pair_types->first << " " << first_pair_types->second << std::endl; // } // Array::const_iterator< Vector > first_pair = pairs.begin(2); // Array::const_iterator< Vector > last_pair = pairs.end(2); // Array & nb_neigh_1 = nb_neighbors(first_pair_types->first, ghost_type1); // Array & nb_neigh_2 = nb_neighbors(first_pair_types->second, ghost_type2); // for(;first_pair != last_pair; ++first_pair) { // UInt q1 = (*first_pair)(0); // UInt q2 = (*first_pair)(1); // ++(nb_neigh_1(q1)); // if(q1 != q2) ++(nb_neigh_2(q2)); // } // } // Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type1); // Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type1); // UInt nb_quads = 0; // Real sum_nb_neig = 0; // UInt max_nb_neig = 0; // UInt min_nb_neig = std::numeric_limits::max(); // for(; it != last_type; ++it) { // Array & nb_neighor = nb_neighbors(*it, ghost_type1); // Array::iterator nb_neigh = nb_neighor.begin(); // Array::iterator end_neigh = nb_neighor.end(); // for (; nb_neigh != end_neigh; ++nb_neigh, ++nb_quads) { // UInt nb = *nb_neigh; // sum_nb_neig += nb; // max_nb_neig = std::max(max_nb_neig, nb); // min_nb_neig = std::min(min_nb_neig, nb); // } // } // comm.allReduce(&nb_quads, 1, _so_sum); // comm.allReduce(&sum_nb_neig, 1, _so_sum); // comm.allReduce(&max_nb_neig, 1, _so_max); // comm.allReduce(&min_nb_neig, 1, _so_min); // if(prank == 0) { // pout << ghost_type2 << " "; // pout << "Nb quadrature points: " << nb_quads << std::endl; // Real mean_nb_neig = sum_nb_neig / Real(nb_quads); // pout << ghost_type2 << " "; // pout << "Average nb neighbors: " << mean_nb_neig << "(" << sum_nb_neig << ")" << std::endl; // pout << ghost_type2 << " "; // pout << "Max nb neighbors: " << max_nb_neig << std::endl; // pout << ghost_type2 << " "; // pout << "Min nb neighbors: " << min_nb_neig << std::endl; // } // } // pout.close(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline UInt MaterialNonLocal::getNbDataForElements(const Array & elements, SynchronizationTag tag) const { UInt nb_quadrature_points = this->getModel().getNbQuadraturePoints(elements); UInt size = 0; if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; size += non_local_variable.nb_component * sizeof(Real) * nb_quadrature_points; } } size += weight_func->getNbDataForElements(elements, tag); return size; } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; this->packElementDataHelper(*non_local_variable.local, buffer, elements); } } weight_func->packElementData(buffer, elements, tag); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) { if(tag == _gst_mnl_for_average) { typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; this->unpackElementDataHelper(*non_local_variable.local, buffer, elements); } } weight_func->unpackElementData(buffer, elements, tag); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::onElementsAdded(const Array & element_list) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ERROR("This is a case not taken into account!!!"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { AKANTU_DEBUG_IN(); // Change the pairs in new global numbering for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::iterator first_pair = pair_list[ghost_type2].begin(); PairList::iterator last_pair = pair_list[ghost_type2].end(); // Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); for(;first_pair != last_pair; ++first_pair) { QuadraturePoint & q1 = first_pair->first; QuadraturePoint gq1 = this->convertToGlobalPoint(q1); q1 = gq1; if(new_numbering.exists(q1.type, q1.ghost_type)) { UInt q1_new_el = new_numbering(q1.type, q1.ghost_type)(gq1.element); AKANTU_DEBUG_ASSERT(q1_new_el != UInt(-1), "A local quadrature_point as been removed instead of just being renumbered"); q1.element = q1_new_el; } QuadraturePoint & q2 = first_pair->second; QuadraturePoint gq2 = this->convertToGlobalPoint(q2); q2 = gq2; if(new_numbering.exists(q2.type, q2.ghost_type)) { UInt q2_new_el = new_numbering(q2.type, q2.ghost_type)(gq2.element); AKANTU_DEBUG_ASSERT(q2_new_el != UInt(-1), "A local quadrature_point as been removed instead of just being renumbered"); q2.element = q2_new_el; } } } // Change the material numbering Material::onElementsRemoved(element_list, new_numbering, event); // Change back the pairs to the new material numbering for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; PairList::iterator first_pair = pair_list[ghost_type2].begin(); PairList::iterator last_pair = pair_list[ghost_type2].end(); // Array::vector_iterator weight_it = pair_weight[ghost_type2]->begin(2); for(;first_pair != last_pair; ++first_pair) { first_pair->first = this->convertToLocalPoint(first_pair->first ); first_pair->second = this->convertToLocalPoint(first_pair->second); } } AKANTU_DEBUG_OUT(); } diff --git a/src/model/solid_mechanics/materials/weight_function.hh b/src/model/solid_mechanics/materials/weight_function.hh index f62661bcf..4fc5fb834 100644 --- a/src/model/solid_mechanics/materials/weight_function.hh +++ b/src/model/solid_mechanics/materials/weight_function.hh @@ -1,373 +1,373 @@ /** * @file weight_function.hh * * @author Nicolas Richart * @author Cyprien Wolff * * @date creation: Fri Apr 13 2012 * @date last modification: Thu Jun 05 2014 * * @brief Weight functions for non local materials * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" #include "solid_mechanics_model.hh" #include "parsable.hh" #include #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #include #endif /* -------------------------------------------------------------------------- */ #include - +#include "material_damage.hh" #ifndef __AKANTU_WEIGHT_FUNCTION_HH__ #define __AKANTU_WEIGHT_FUNCTION_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Normal weight function */ /* -------------------------------------------------------------------------- */ template class BaseWeightFunction : public Parsable { public: BaseWeightFunction(Material & material, const std::string & type = "base") : Parsable(_st_non_local, "weight_function:" + type), material(material), type(type) { this->registerParam("radius" , R , 100., _pat_parsable | _pat_readable , "Non local radius"); this->registerParam("update_rate" , update_rate, 0U , _pat_parsmod, "Update frequency"); } virtual ~BaseWeightFunction() {} virtual void init() { R2 = R * R; }; virtual void updateInternals(__attribute__((unused)) const ElementTypeMapArray & quadrature_points_coordinates) {}; /* ------------------------------------------------------------------------ */ inline void setRadius(Real radius) { R = radius; R2 = R * R; } /* ------------------------------------------------------------------------ */ inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, __attribute__((unused)) ElementType type2, __attribute__((unused)) GhostType ghost_type2) { } /* ------------------------------------------------------------------------ */ inline Real operator()(Real r, const __attribute__((unused)) QuadraturePoint & q1, const __attribute__((unused)) QuadraturePoint & q2) { Real w = 0; if(r <= R) { Real alpha = (1. - r*r / R2); w = alpha * alpha; // *weight = 1 - sqrt(r / radius); } return w; } void printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "WeightFunction " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } public: Real getRadius() { return R; } UInt getUpdateRate() { return update_rate; } public: virtual UInt getNbDataForElements(__attribute__((unused)) const Array & elements, __attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline void packElementData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Array & elements, __attribute__((unused)) SynchronizationTag tag) const {} virtual inline void unpackElementData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Array & elements, __attribute__((unused)) SynchronizationTag tag) {} protected: Material & material; Real R; Real R2; UInt update_rate; const std::string type; }; /* -------------------------------------------------------------------------- */ /* Damage weight function */ /* -------------------------------------------------------------------------- */ template class DamagedWeightFunction : public BaseWeightFunction { public: DamagedWeightFunction(Material & material) : BaseWeightFunction(material, "damaged") { - AKANTU_DEBUG_ASSERT(dynamic_cast *>(&material) != NULL, "This weight function works only with damage materials!"); + //AKANTU_DEBUG_ASSERT(dynamic_cast *>(&material) != NULL, "This weight function works only with damage materials!"); } inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_damage = &dynamic_cast &>(this->material).getDamage(type2, ghost_type2); } inline Real operator()(Real r, const __attribute__((unused)) QuadraturePoint & q1, const QuadraturePoint & q2) { UInt quad = q2.global_num; Real D = (*selected_damage)(quad); Real Radius_t = 0; Real Radius_init = this->R2; // if(D <= 0.5) // { // Radius_t = 2*D*Radius_init; // } // else // { // Radius_t = 2*Radius_init*(1-D); // } // Radius_t = Radius_init*(1-D); Radius_init *= Radius_init; Radius_t *= Radius_t; if(Radius_t < Math::getTolerance()) { Radius_t = 0.001*Radius_init; } Real expb = (2*std::log(0.51))/(std::log(1.0-0.49*Radius_t/Radius_init)); Int expb_floor=std::floor(expb); Real b = expb_floor + expb_floor%2; Real alpha = std::max(0., 1. - r*r / Radius_init); Real w = std::pow(alpha,b); return w; } private: const Array * selected_damage; }; /* -------------------------------------------------------------------------- */ /* Remove damaged weight function */ /* -------------------------------------------------------------------------- */ template class RemoveDamagedWeightFunction : public BaseWeightFunction { public: RemoveDamagedWeightFunction(Material & material) : BaseWeightFunction(material, "remove_damaged") { AKANTU_DEBUG_ASSERT(dynamic_cast *>(&material) != NULL, "This weight function works only with damage materials!"); this->registerParam("damage_limit", this->damage_limit, 1., _pat_parsable, "Damage Threshold"); } inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { MaterialDamage & mat = dynamic_cast &>(this->material); selected_damage = &mat.getDamage(type2, ghost_type2); } inline Real operator()(Real r, const __attribute__((unused)) QuadraturePoint & q1, const QuadraturePoint & q2) { UInt quad = q2.global_num; if(q1 == q2) return 1.; Real D = (*selected_damage)(quad); Real w = 0.; if(D < damage_limit) { Real alpha = std::max(0., 1. - r*r / this->R2); w = alpha * alpha; } return w; } virtual UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_weight) return this->material.getModel().getNbQuadraturePoints(elements) * sizeof(Real); return 0; } virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_weight) { ElementTypeMapArray & damage = this->material.getInternal("damage"); this->material.packElementDataHelper(damage, buffer, elements); } } virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) { if(tag == _gst_mnl_weight) { ElementTypeMapArray & damage = this->material.getInternal("damage"); this->material.unpackElementDataHelper(damage, buffer, elements); } } private: /// limit at which a point is considered as complitely broken Real damage_limit; /// internal pointer to the current damage vector const Array * selected_damage; }; /* -------------------------------------------------------------------------- */ /* Remove damaged with damage rate weight function */ /* -------------------------------------------------------------------------- */ template class RemoveDamagedWithDamageRateWeightFunction : public BaseWeightFunction { public: RemoveDamagedWithDamageRateWeightFunction(Material & material) : BaseWeightFunction(material, "remove_damage_with_damage_rate") { this->registerParam("damage_limit", this->damage_limit_with_damage_rate, 1, _pat_parsable, "Damage Threshold"); } inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_damage_with_damage_rate = &(this->material.getArray("damage",type2, ghost_type2)); selected_damage_rate_with_damage_rate = &(this->material.getArray("damage-rate",type2, ghost_type2)); } inline Real operator()(Real r, const __attribute__((unused)) QuadraturePoint & q1, const QuadraturePoint & q2) { UInt quad = q2.global_num; if(q1.global_num == quad) return 1.; Real D = (*selected_damage_with_damage_rate)(quad); Real w = 0.; Real alphaexp = 1.; Real betaexp = 2.; if(D < damage_limit_with_damage_rate) { Real alpha = std::max(0., 1. - pow((r*r / this->R2),alphaexp)); w = pow(alpha, betaexp); } return w; } private: /// limit at which a point is considered as complitely broken Real damage_limit_with_damage_rate; /// internal pointer to the current damage vector const Array * selected_damage_with_damage_rate; /// internal pointer to the current damage rate vector const Array * selected_damage_rate_with_damage_rate; }; /* -------------------------------------------------------------------------- */ /* Stress Based Weight */ /* -------------------------------------------------------------------------- */ template class StressBasedWeightFunction : public BaseWeightFunction { public: StressBasedWeightFunction(Material & material); void init(); virtual void updateInternals(__attribute__((unused)) const ElementTypeMapArray & quadrature_points_coordinates) { updatePrincipalStress(_not_ghost); updatePrincipalStress(_ghost); }; void updatePrincipalStress(GhostType ghost_type); inline void updateQuadraturePointsCoordinates(ElementTypeMapArray & quadrature_points_coordinates); inline void selectType(ElementType type1, GhostType ghost_type1, ElementType type2, GhostType ghost_type2); inline Real operator()(Real r, const QuadraturePoint & q1, const QuadraturePoint & q2); inline Real computeRhoSquare(Real r, Vector & eigs, Matrix & eigenvects, Vector & x_s); private: Real ft; InternalField stress_diag; Array * selected_stress_diag; InternalField stress_base; Array * selected_stress_base; // InternalField quadrature_points_coordinates; Array * selected_position_1; Array * selected_position_2; InternalField characteristic_size; Array * selected_characteristic_size; }; template inline std::ostream & operator <<(std::ostream & stream, const BaseWeightFunction & _this) { _this.printself(stream); return stream; } #include "weight_function_tmpl.hh" __END_AKANTU__ #endif /* __AKANTU_WEIGHT_FUNCTION_HH__ */ diff --git a/src/solver/petsc_matrix.cc b/src/solver/petsc_matrix.cc index 35b40980f..6204a0c13 100644 --- a/src/solver/petsc_matrix.cc +++ b/src/solver/petsc_matrix.cc @@ -1,630 +1,636 @@ /** * @file petsc_matrix.cc * @author Aurelia Cuba Ramos * @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 . * */ /* -------------------------------------------------------------------------- */ #include "petsc_matrix.hh" #include "static_communicator.hh" #include "static_communicator_mpi.hh" #include "mpi_type_wrapper.hh" #include "dof_synchronizer.hh" #include "petsc_wrapper.hh" /* -------------------------------------------------------------------------- */ #include #include __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, const ID & id, const MemoryID & memory_id) : SparseMatrix(size, sparse_matrix_type, id, memory_id), petsc_matrix_wrapper(new PETScMatrixWrapper), d_nnz(0,1,"dnnz"), o_nnz(0,1,"onnz"), first_global_index(0), is_petsc_matrix_initialized(false) { AKANTU_DEBUG_IN(); StaticSolver::getStaticSolver().registerEventHandler(*this); this->offset = 0; this->init(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PETScMatrix::PETScMatrix(const SparseMatrix & matrix, const ID & id, const MemoryID & memory_id) : SparseMatrix(matrix, id, memory_id), petsc_matrix_wrapper(new PETScMatrixWrapper), d_nnz(0,1,"dnnz"), o_nnz(0,1,"onnz"), first_global_index(0), is_petsc_matrix_initialized(false) { // AKANTU_DEBUG_IN(); // StaticSolver::getStaticSolver().registerEventHandler(*this); // this->offset = 0; // this->init(); // AKANTU_DEBUG_OUT(); AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ PETScMatrix::~PETScMatrix() { AKANTU_DEBUG_IN(); /// destroy all the PETSc data structures used for this matrix this->destroyInternalData(); + StaticSolver::getStaticSolver().unregisterEventHandler(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PETScMatrix::init() { AKANTU_DEBUG_IN(); /// set the communicator that should be used by PETSc #if defined(AKANTU_USE_MPI) StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); const StaticCommunicatorMPI & mpi_st_comm = dynamic_cast(comm.getRealStaticCommunicator()); this->petsc_matrix_wrapper->communicator = mpi_st_comm.getMPITypeWrapper().getMPICommunicator(); #else this->petsc_matrix_wrapper->communicator = PETSC_COMM_SELF; #endif PetscErrorCode ierr; /// create the PETSc matrix object ierr = MatCreate(this->petsc_matrix_wrapper->communicator, &(this->petsc_matrix_wrapper->mat)); CHKERRXX(ierr); /** * Set the matrix type * @todo PETSc does currently not support a straightforward way to * apply Dirichlet boundary conditions for MPISBAIJ * matrices. Therefore always the entire matrix is allocated. It * would be possible to use MATSBAIJ for sequential matrices in case * memory becomes critical. Also, block matrices would give a much * better performance. Modify this in the future! */ ierr = MatSetType(this->petsc_matrix_wrapper->mat, MATAIJ); CHKERRXX(ierr); this->is_petsc_matrix_initialized = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * With this method each processor computes the dimensions of the * local matrix, i.e. the part of the global matrix it is storing. * @param dof_synchronizer dof synchronizer that maps the local * dofs to the global dofs and the equation numbers, i.e., the * position at which the dof is assembled in the matrix */ void PETScMatrix::setSize() { AKANTU_DEBUG_IN(); PetscErrorCode ierr; /// find the number of dofs corresponding to master or local nodes UInt nb_dofs = this->dof_synchronizer->getNbDOFs(); UInt nb_local_master_dofs = 0; /// create array to store the global equation number of all local and master dofs Array local_master_eq_nbs(nb_dofs); Array::scalar_iterator it_eq_nb = local_master_eq_nbs.begin(); /// get the pointer to the global equation number array Int * eq_nb_val = this->dof_synchronizer->getGlobalDOFEquationNumbers().storage(); for (UInt i = 0; i dof_synchronizer->isLocalOrMasterDOF(i) ) { *it_eq_nb = eq_nb_val[i]; ++it_eq_nb; ++nb_local_master_dofs; } } local_master_eq_nbs.resize(nb_local_master_dofs); /// set the local size this->local_size = nb_local_master_dofs; /// resize PETSc matrix #if defined(AKANTU_USE_MPI) ierr = MatSetSizes(this->petsc_matrix_wrapper->mat, this->local_size, this->local_size, this->size, this->size); CHKERRXX(ierr); #else ierr = MatSetSizes(this->petsc_matrix_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(); } /* -------------------------------------------------------------------------- */ /** * This method generates a mapping from the global Akantu equation * numbering to the global PETSc dof ordering * @param local_master_eq_nbs_ptr Int pointer to the array of equation * numbers of all local or master dofs, i.e. the row indices of the * local matrix */ void PETScMatrix::createGlobalAkantuToPETScMap(Int* local_master_eq_nbs_ptr) { AKANTU_DEBUG_IN(); PetscErrorCode ierr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt rank = comm.whoAmI(); //initialize vector to store the number of local and master nodes that are assigned to each processor Vector master_local_ndofs_per_proc(nb_proc); /// store the nb of master and local dofs on each processor master_local_ndofs_per_proc(rank) = this->local_size; /// 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 /// determine the PETSc-index for the first dof on each processor for (UInt i = 0; i < rank; ++i) { this->first_global_index += master_local_ndofs_per_proc(i); } /// create array for petsc ordering Array petsc_dofs(this->local_size); Array::scalar_iterator it_petsc_dofs = petsc_dofs.begin(); for (Int i = this->first_global_index; i < this->first_global_index + this->local_size; ++i, ++it_petsc_dofs) { *it_petsc_dofs = i; } ierr = AOCreateBasic(this->petsc_matrix_wrapper->communicator, this->local_size, local_master_eq_nbs_ptr, petsc_dofs.storage(), &(this->petsc_matrix_wrapper->ao)); CHKERRXX(ierr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void PETScMatrix::createLocalAkantuToPETScMap(const DOFSynchronizer & dof_synchronizer) { // AKANTU_DEBUG_IN(); // AKANTU_DEBUG_ASSERT(this->petsc_matrix_wrapper->ao != NULL, // "You should first create a mapping from the global" // << " Akantu numbering to the global PETSc numbering"); // PetscErrorCode ierr; // this->dof_synchronizer = &const_cast(dof_synchronizer); // /// get the number of dofs // Int nb_dofs = dof_synchronizer.getNbDOFs(); // /// get the global equation numbers for each node // Array global_dof_equation_numbers = dof_synchronizer.getGlobalDOFEquationNumbers(); // /// map the global dof equation numbers to the corresponding PETSc ordering // ierr = AOApplicationToPETSc(this->petsc_matrix_wrapper->ao, nb_dofs, // global_dof_equation_numbers.storage()); CHKERRXX(ierr); // /// create the mapping from the local Akantu ordering to the global PETSc ordering // ierr = ISLocalToGlobalMappingCreate(this->petsc_matrix_wrapper->communicator, // 1, nb_dofs, global_dof_equation_numbers.storage(), // PETSC_COPY_VALUES, &(this->petsc_matrix_wrapper-mapping)); CHKERRXX(ierr); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ /** * This function creates the non-zero pattern of the PETSc matrix. In * PETSc the parallel matrix is partitioned across processors such * that the first m0 rows belong to process 0, the next m1 rows belong * to process 1, the next m2 rows belong to process 2 etc.. where * m0,m1,m2,.. are the input parameter 'm'. i.e each processor stores * values corresponding to [m x N] submatrix * (http://www.mcs.anl.gov/petsc/). * @param mesh mesh discretizing the domain we want to analyze * @param dof_synchronizer dof synchronizer that maps the local * dofs to the global dofs and the equation numbers, i.e., the * position at which the dof is assembled in the matrix */ void PETScMatrix::buildProfile(const Mesh & mesh, const DOFSynchronizer & dof_synchronizer, UInt nb_degree_of_freedom) { AKANTU_DEBUG_IN(); //clearProfile(); this->dof_synchronizer = &const_cast(dof_synchronizer); this->setSize(); PetscErrorCode ierr; /// 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); /// set arrays to zero everywhere this->d_nnz.set(0); this->o_nnz.set(0); // if(irn_jcn_to_k) delete irn_jcn_to_k; // irn_jcn_to_k = new std::map, UInt>; coordinate_list_map::iterator irn_jcn_k_it; Int * eq_nb_val = dof_synchronizer.getGlobalDOFEquationNumbers().storage(); UInt nb_global_dofs = dof_synchronizer.getNbGlobalDOFs(); /// 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; 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); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, ghost_type); 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, ghost_type).storage(); Int * local_eq_nb_val = new Int[nb_degree_of_freedom * nb_nodes_per_element]; 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) { /** * !!!!!!Careful!!!!!! This is a ugly fix. @todo this is a * very ugly fix, because the offset for the global * equation number, where the dof will be assembled, is * hardcoded. In the future a class dof manager has to be * added to Akantu to handle the mapping between the dofs * and the equation numbers * */ *tmp_local_eq_nb_val++ = eq_nb_val[n * nb_degree_of_freedom + d] - (dof_synchronizer.isPureGhostDOF(n * nb_degree_of_freedom + d) ? nb_global_dofs : 0); } } for (UInt i = 0; i < size_mat; ++i) { Int c_irn = local_eq_nb_val[i]; UInt j_start = 0; for (UInt j = j_start; j < size_mat; ++j) { Int c_jcn = local_eq_nb_val[j]; Array index_pair(2,1); index_pair(0) = c_irn; index_pair(1) = c_jcn; AOApplicationToPetsc(this->petsc_matrix_wrapper->ao, 2, index_pair.storage()); if (index_pair(0) >= first_global_index && index_pair(0) < first_global_index + this->local_size) { KeyCOO irn_jcn = keyPETSc(c_irn, c_jcn); 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; // check if node is slave node if (index_pair(1) >= first_global_index && index_pair(1) < first_global_index + this->local_size) this->d_nnz(index_pair(0) - first_global_index) += 1; else this->o_nnz(index_pair(0) - first_global_index) += 1; nb_non_zero++; } } } } conn_val += nb_nodes_per_element; } delete [] local_eq_nb_val; } } // /// 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++; // } // } // } // 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_matrix_wrapper->mat, &type); - std::cout << "the matrix type is: " << type << std::endl; + /// std::cout << "the matrix type is: " << type << std::endl; /** * PETSc will use only one of the following preallocation commands * depending on the matrix type and ignore the rest. Note that for * the block matrix format a block size of 1 is used. This might * result in bad performance. @todo For better results implement * buildProfile() with larger block size. * */ /// build profile: if (strcmp(type, MATSEQAIJ) == 0) { ierr = MatSeqAIJSetPreallocation(this->petsc_matrix_wrapper->mat, 0, d_nnz.storage()); CHKERRXX(ierr); } else if ((strcmp(type, MATMPIAIJ) == 0)) { ierr = MatMPIAIJSetPreallocation(this->petsc_matrix_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"); } //ierr = MatSeqSBAIJSetPreallocation(this->petsc_matrix_wrapper->mat, 1, // 0, d_nnz.storage()); CHKERRXX(ierr); if (this->sparse_matrix_type==_symmetric) { /// set flag for symmetry to enable ICC/Cholesky preconditioner ierr = MatSetOption(this->petsc_matrix_wrapper->mat, MAT_SYMMETRIC, PETSC_TRUE); CHKERRXX(ierr); /// set flag for symmetric positive definite ierr = MatSetOption(this->petsc_matrix_wrapper->mat, MAT_SPD, PETSC_TRUE); CHKERRXX(ierr); } /// once the profile has been build ignore any new nonzero locations ierr = MatSetOption(this->petsc_matrix_wrapper->mat, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE); CHKERRXX(ierr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Method to save the nonzero pattern and the values stored at each position * @param filename name of the file in which the information will be stored */ void PETScMatrix::saveMatrix(const std::string & filename) const{ AKANTU_DEBUG_IN(); PetscErrorCode ierr; /// create Petsc viewer PetscViewer viewer; ierr = PetscViewerASCIIOpen(this->petsc_matrix_wrapper->communicator, filename.c_str(), &viewer); CHKERRXX(ierr); /// set the format PetscViewerSetFormat(viewer, PETSC_VIEWER_DEFAULT); CHKERRXX(ierr); /// save the matrix + /// @todo Write should be done in serial -> might cause problems ierr = MatView(this->petsc_matrix_wrapper->mat, viewer); CHKERRXX(ierr); /// destroy the viewer ierr = PetscViewerDestroy(&viewer); CHKERRXX(ierr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Method to add an Akantu sparse matrix to the PETSc matrix * @param matrix Akantu sparse matrix to be added * @param alpha the factor specifying how many times the matrix should be added */ void PETScMatrix::add(const SparseMatrix & matrix, Real alpha) { PetscErrorCode ierr; // AKANTU_DEBUG_ASSERT(nb_non_zero == matrix.getNbNonZero(), // "The two matrix don't have the same profiles"); + Real val_to_add = 0; for (UInt n = 0; n < matrix.getNbNonZero(); ++n) { Array index(2,1); - index(0) = matrix.getIRN()(n)-matrix.getOffset(); - index(1) = matrix.getJCN()(n)-matrix.getOffset(); + UInt mat_to_add_offset = matrix.getOffset(); + index(0) = matrix.getIRN()(n)-mat_to_add_offset; + index(1) = matrix.getJCN()(n)-mat_to_add_offset; AOApplicationToPetsc(this->petsc_matrix_wrapper->ao, 2, index.storage()); if (this->sparse_matrix_type == _symmetric && index(0) > index(1)) std::swap(index(0), index(1)); + + val_to_add = matrix.getA()(n) * alpha; /// MatSetValue might be very slow for MATBAIJ, might need to use MatSetValuesBlocked - ierr = MatSetValue(this->petsc_matrix_wrapper->mat, index(0), index(1), matrix.getA()(n) * alpha, ADD_VALUES); CHKERRXX(ierr); + ierr = MatSetValue(this->petsc_matrix_wrapper->mat, index(0), index(1), val_to_add, ADD_VALUES); CHKERRXX(ierr); /// chek if sparse matrix to be added is symmetric. In this case /// the value also needs to be added at the transposed location in /// the matrix because PETSc is using the full profile, also for symmetric matrices if (matrix.getSparseMatrixType() == _symmetric && index(0) != index(1)) - ierr = MatSetValue(this->petsc_matrix_wrapper->mat, index(1), index(0), matrix.getA()(n) * alpha, ADD_VALUES); CHKERRXX(ierr); + ierr = MatSetValue(this->petsc_matrix_wrapper->mat, index(1), index(0), val_to_add, ADD_VALUES); CHKERRXX(ierr); } this->performAssembly(); } /* -------------------------------------------------------------------------- */ /** * Method to add another PETSc matrix to this PETSc matrix * @param matrix PETSc matrix to be added * @param alpha the factor specifying how many times the matrix should be added */ void PETScMatrix::add(const PETScMatrix & matrix, Real alpha) { PetscErrorCode ierr; ierr = MatAXPY(this->petsc_matrix_wrapper->mat, alpha, matrix.petsc_matrix_wrapper->mat, SAME_NONZERO_PATTERN); CHKERRXX(ierr); this->performAssembly(); } /* -------------------------------------------------------------------------- */ /** * MatSetValues() generally caches the values. The matrix is ready to * use only after MatAssemblyBegin() and MatAssemblyEnd() have been * called. (http://www.mcs.anl.gov/petsc/) */ void PETScMatrix::performAssembly() { PetscErrorCode ierr; ierr = MatAssemblyBegin(this->petsc_matrix_wrapper->mat, MAT_FINAL_ASSEMBLY); CHKERRXX(ierr); ierr = MatAssemblyEnd(this->petsc_matrix_wrapper->mat, MAT_FINAL_ASSEMBLY); CHKERRXX(ierr); } /* -------------------------------------------------------------------------- */ /** * Method is called when the static solver is destroyed, just before * PETSc is finalized. So all PETSc objects need to be destroyed at * this point. */ void PETScMatrix::beforeStaticSolverDestroy() { AKANTU_DEBUG_IN(); try{ this->destroyInternalData(); } catch(...) {} AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// destroy all the PETSc data structures used for this matrix void PETScMatrix::destroyInternalData() { AKANTU_DEBUG_IN(); if(this->is_petsc_matrix_initialized) { PetscErrorCode ierr; ierr = MatDestroy(&(this->petsc_matrix_wrapper->mat)); CHKERRXX(ierr); delete petsc_matrix_wrapper; this->is_petsc_matrix_initialized = false; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// access K(i, j). Works only for dofs on this processor!!!! Real PETScMatrix::operator()(UInt i, UInt j) const{ AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(this->dof_synchronizer->isLocalOrMasterDOF(i) && this->dof_synchronizer->isLocalOrMasterDOF(j), "Operator works only for dofs on this processor"); Array index(2,1); index(0) = this->dof_synchronizer->getDOFGlobalID(i); index(1) = this->dof_synchronizer->getDOFGlobalID(j); AOApplicationToPetsc(this->petsc_matrix_wrapper->ao, 2, index.storage()); Real value = 0; PetscErrorCode ierr; /// @todo MatGetValue might be very slow for MATBAIJ, might need to use MatGetValuesBlocked ierr = MatGetValues(this->petsc_matrix_wrapper->mat, 1, &index(0), 1, &index(1), &value); CHKERRXX(ierr); AKANTU_DEBUG_OUT(); return value; } /* -------------------------------------------------------------------------- */ /** * Apply Dirichlet boundary conditions by zeroing the rows and columns which correspond to blocked dofs * @param boundary array of booleans which are true if the dof at this position is blocked * @param block_val the value in the diagonal entry of blocked rows */ void PETScMatrix::applyBoundary(const Array & boundary, Real block_val) { AKANTU_DEBUG_IN(); PetscErrorCode ierr; /// get the global equation numbers to find the rows that need to be zeroed for the blocked dofs Int * eq_nb_val = dof_synchronizer->getGlobalDOFEquationNumbers().storage(); /// every processor calls the MatSetZero() only for his local or master dofs. This assures that not two processors or more try to zero the same row - Int nb_blocked_local_master_eq_nb = 0; - Array blocked_local_master_eq_nb(this->local_size / boundary.getNbComponent(), boundary.getNbComponent()); - Int * blocked_lm_eq_nb_ptr = blocked_local_master_eq_nb.storage(); UInt nb_component = boundary.getNbComponent(); UInt size = boundary.getSize(); + Int nb_blocked_local_master_eq_nb = 0; + Array blocked_local_master_eq_nb(this->local_size / nb_component, nb_component); + Int * blocked_lm_eq_nb_ptr = blocked_local_master_eq_nb.storage(); for (UInt i = 0; i < size; ++i) { for (UInt j = 0; j < nb_component; ++j) { UInt local_dof = i * nb_component + j; if (boundary(i, j) == true && this->dof_synchronizer->isLocalOrMasterDOF(local_dof)) { Int global_eq_nb = *eq_nb_val; *blocked_lm_eq_nb_ptr = global_eq_nb; ++nb_blocked_local_master_eq_nb; ++blocked_lm_eq_nb_ptr; } ++eq_nb_val; } } - blocked_local_master_eq_nb.resize(nb_blocked_local_master_eq_nb/boundary.getNbComponent()); + blocked_local_master_eq_nb.resize(nb_blocked_local_master_eq_nb/nb_component); ierr = AOApplicationToPetsc(this->petsc_matrix_wrapper->ao, nb_blocked_local_master_eq_nb, blocked_local_master_eq_nb.storage() ); CHKERRXX(ierr); ierr = MatZeroRowsColumns(this->petsc_matrix_wrapper->mat, nb_blocked_local_master_eq_nb, blocked_local_master_eq_nb.storage(), block_val, 0, 0); CHKERRXX(ierr); this->performAssembly(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// set all entries to zero while keeping the same nonzero pattern void PETScMatrix::clear() { MatZeroEntries(this->petsc_matrix_wrapper->mat); } __END_AKANTU__ diff --git a/src/solver/solver_petsc.cc b/src/solver/solver_petsc.cc index dc9f44e7b..635c10ae6 100644 --- a/src/solver/solver_petsc.cc +++ b/src/solver/solver_petsc.cc @@ -1,1106 +1,1109 @@ /** * @file solver_petsc.cc * * @author Nicolas Richart # @author Alejandro M. Aragón * @author Aurelia Cuba Ramos * * @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 . * */ /* -------------------------------------------------------------------------- */ #include #include "solver_petsc.hh" #include "petsc_wrapper.hh" #include "petsc_matrix.hh" #include "static_communicator.hh" #include "static_communicator_mpi.hh" #include "mpi_type_wrapper.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolverPETSc::SolverPETSc(SparseMatrix & matrix, const ID & id, const MemoryID & memory_id) : Solver(matrix, id, memory_id), is_petsc_data_initialized(false), petsc_solver_wrapper(new PETScSolverWrapper) { } /* -------------------------------------------------------------------------- */ SolverPETSc::~SolverPETSc() { AKANTU_DEBUG_IN(); this->destroyInternalData(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverPETSc::destroyInternalData() { AKANTU_DEBUG_IN(); if(this->is_petsc_data_initialized) { PetscErrorCode ierr; ierr = KSPDestroy(&(this->petsc_solver_wrapper->ksp)); CHKERRXX(ierr); ierr = VecDestroy(&(this->petsc_solver_wrapper->rhs)); CHKERRXX(ierr); ierr = VecDestroy(&(this->petsc_solver_wrapper->solution)); CHKERRXX(ierr); delete petsc_solver_wrapper; this->is_petsc_data_initialized = false; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverPETSc::initialize(SolverOptions & options) { AKANTU_DEBUG_IN(); #if defined(AKANTU_USE_MPI) StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); const StaticCommunicatorMPI & mpi_st_comm = dynamic_cast(comm.getRealStaticCommunicator()); this->petsc_solver_wrapper->communicator = mpi_st_comm.getMPITypeWrapper().getMPICommunicator(); #else this->petsc_solver_wrapper->communicator = PETSC_COMM_SELF; #endif PetscErrorCode ierr; PETScMatrix * petsc_matrix = static_cast(this->matrix); /// create a solver context ierr = KSPCreate(this->petsc_solver_wrapper->communicator, &(this->petsc_solver_wrapper->ksp)); CHKERRXX(ierr); + /// create the PETSc vector for the right-hand side ierr = VecCreate(this->petsc_solver_wrapper->communicator, &(this->petsc_solver_wrapper->rhs)); CHKERRXX(ierr); ierr = VecSetSizes((this->petsc_solver_wrapper->rhs), petsc_matrix->getLocalSize(), petsc_matrix->getSize()); CHKERRXX(ierr); ierr = VecSetFromOptions((this->petsc_solver_wrapper->rhs)); CHKERRXX(ierr); + /// create the PETSc vector for the solution ierr = VecDuplicate((this->petsc_solver_wrapper->rhs), &(this->petsc_solver_wrapper->solution)); CHKERRXX(ierr); + /// set the solution to zero ierr = VecZeroEntries(this->petsc_solver_wrapper->solution); this->is_petsc_data_initialized = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverPETSc::setRHS(Array & rhs) { PetscErrorCode ierr; PETScMatrix * petsc_matrix = static_cast(this->matrix); UInt nb_component = rhs.getNbComponent(); UInt size = rhs.getSize(); for (UInt i = 0; i < size; ++i) { for (UInt j = 0; j < nb_component; ++j) { UInt i_local = i * nb_component + j; if (this->matrix->getDOFSynchronizer().isLocalOrMasterDOF(i_local)) { Int i_global = this->matrix->getDOFSynchronizer().getDOFGlobalID(i_local); AOApplicationToPetsc(petsc_matrix->getPETScMatrixWrapper()->ao, 1, &(i_global) ); ierr = VecSetValue((this->petsc_solver_wrapper->rhs), i_global, rhs(i,j), INSERT_VALUES); CHKERRXX(ierr); } } } ierr = VecAssemblyBegin(this->petsc_solver_wrapper->rhs); CHKERRXX(ierr); ierr = VecAssemblyEnd(this->petsc_solver_wrapper->rhs); CHKERRXX(ierr); /// ierr = VecCopy((this->petsc_solver_wrapper->rhs), (this->petsc_solver_wrapper->solution)); CHKERRXX(ierr); } /* -------------------------------------------------------------------------- */ void SolverPETSc::solve() { AKANTU_DEBUG_IN(); PetscErrorCode ierr; ierr = KSPSolve(this->petsc_solver_wrapper->ksp, this->petsc_solver_wrapper->rhs, this->petsc_solver_wrapper->solution); CHKERRXX(ierr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverPETSc::solve(Array & solution) { AKANTU_DEBUG_IN(); this->solution = &solution; this->solve(); PetscErrorCode ierr; PETScMatrix * petsc_matrix = static_cast(this->matrix); // ierr = VecGetOwnershipRange(this->petsc_solver_wrapper->solution, solution_begin, solution_end); CHKERRXX(ierr); // ierr = VecGetArray(this->petsc_solver_wrapper->solution, PetscScalar **array); CHKERRXX(ierr); UInt nb_component = solution.getNbComponent(); UInt size = solution.getSize(); for (UInt i = 0; i < size; ++i) { for (UInt j = 0; j < nb_component; ++j) { UInt i_local = i * nb_component + j; if (this->matrix->getDOFSynchronizer().isLocalOrMasterDOF(i_local)) { Int i_global = this->matrix->getDOFSynchronizer().getDOFGlobalID(i_local); ierr = AOApplicationToPetsc(petsc_matrix->getPETScMatrixWrapper()->ao, 1, &(i_global) ); CHKERRXX(ierr); ierr = VecGetValues(this->petsc_solver_wrapper->solution, 1, &i_global, &solution(i,j)); CHKERRXX(ierr); } } } synch_registry->synchronize(_gst_solver_solution); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverPETSc::setOperators() { AKANTU_DEBUG_IN(); PetscErrorCode ierr; /// set the matrix that defines the linear system and the matrix for preconditioning (here they are the same) PETScMatrix * petsc_matrix = static_cast(this->matrix); #if PETSC_VERSION_MAJOR >= 3 && PETSC_VERSION_MINOR >= 5 ierr = KSPSetOperators(this->petsc_solver_wrapper->ksp, petsc_matrix->getPETScMatrixWrapper()->mat, petsc_matrix->getPETScMatrixWrapper()->mat); CHKERRXX(ierr); #else ierr = KSPSetOperators(this->petsc_solver_wrapper->ksp, petsc_matrix->getPETScMatrixWrapper()->mat, petsc_matrix->getPETScMatrixWrapper()->mat, SAME_NONZERO_PATTERN); CHKERRXX(ierr); #endif /// If this is not called the solution vector is zeroed in the call to KSPSolve(). ierr = KSPSetInitialGuessNonzero(this->petsc_solver_wrapper->ksp, PETSC_TRUE); KSPSetFromOptions(this->petsc_solver_wrapper->ksp); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void finalize_petsc() { // static bool finalized = false; // if (!finalized) { // cout<<"*** INFO *** PETSc is finalizing..."<first; // ierr = VecSetValues(b, 1, &row, &it->second, ADD_VALUES); // } // #ifdef CPPUTILS_VERBOSE // out<<" Right hand side set... "< 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... "< -pc_type -ksp_monitor -ksp_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... "<first; // ierr = VecGetValues(x_all, 1, &row, &val); // if (val != zero) // xx[row] = val; // } // #ifdef CPPUTILS_VERBOSE // out<<" Solution vector copied... "< Kpf_nz = KKpf.non_zero_off_diagonal(); // std::pair 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); // for (sparse_matrix_type::const_hash_iterator it = KKpf.map_.begin(); it != KKpf.map_.end(); ++it) { // // get subscripts // std::pair 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 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... "<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); // // add Kpf*Uf // ierr = MatMult(Kpf, x_, Pf); // // add Kpp*Up // ierr = MatMultAdd(Kpp, Up, Pf, Pp); // #ifdef CPPUTILS_VERBOSE // out<<" Matrices multiplied... "<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... "< 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..."<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..."<(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(&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(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(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(&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 & 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 & 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__ diff --git a/test/test_solver/1D_bar.geo b/test/test_solver/1D_bar.geo new file mode 100644 index 000000000..235f16da2 --- /dev/null +++ b/test/test_solver/1D_bar.geo @@ -0,0 +1,4 @@ +Point(1) = {0, 0, 0, 1}; +Point(2) = {10, 0, 0, 1}; +Line(1) = {1, 2}; +Physical Line(2) = {1}; diff --git a/test/test_solver/CMakeLists.txt b/test/test_solver/CMakeLists.txt index fd1de042c..61ee1e556 100644 --- a/test/test_solver/CMakeLists.txt +++ b/test/test_solver/CMakeLists.txt @@ -1,70 +1,84 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author Nicolas Richart # # @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 . # # @section DESCRIPTION # #=============================================================================== add_mesh(test_solver_mesh triangle.geo 2 1) +add_mesh(test_matrix_mesh square.geo 2 1) +add_mesh(test_solver_petsc_mesh 1D_bar.geo 1 1) register_test(test_sparse_matrix_profile SOURCES test_sparse_matrix_profile.cc DEPENDENCIES test_solver_mesh PACKAGE implicit ) register_test(test_sparse_matrix_assemble SOURCES test_sparse_matrix_assemble.cc DEPENDENCIES test_solver_mesh PACKAGE implicit ) register_test(test_sparse_matrix_product SOURCES test_sparse_matrix_product.cc FILES_TO_COPY bar.msh PACKAGE implicit ) register_test(test_petsc_matrix_profile SOURCES test_petsc_matrix_profile.cc + DEPENDENCIES test_matrix_mesh + PACKAGE petsc + ) + +register_test(test_petsc_matrix_profile_parallel + SOURCES test_petsc_matrix_profile_parallel.cc + DEPENDENCIES test_matrix_mesh + PACKAGE petsc + ) + +register_test(test_petsc_matrix_diagonal + SOURCES test_petsc_matrix_diagonal.cc DEPENDENCIES test_solver_mesh PACKAGE petsc ) register_test(test_petsc_matrix_apply_boundary SOURCES test_petsc_matrix_apply_boundary.cc DEPENDENCIES test_solver_mesh PACKAGE petsc ) register_test(test_solver_petsc SOURCES test_solver_petsc.cc - DEPENDENCIES profile.mtx + DEPENDENCIES test_solver_petsc_mesh PACKAGE petsc ) \ No newline at end of file diff --git a/test/test_solver/square.geo b/test/test_solver/square.geo new file mode 100644 index 000000000..67df205f8 --- /dev/null +++ b/test/test_solver/square.geo @@ -0,0 +1,28 @@ +// Mesh size +h = 1.0; // Top cube + +// Dimensions of top cube +Lx = 1; +Ly = 1; + +// ------------------------------------------ +// Geometry +// ------------------------------------------ + +// Base Cube +Point(101) = { 0.0, 0.0, 0.0, h}; // Bottom Face +Point(102) = { Lx, 0.0, 0.0, h}; // Bottom Face +Point(103) = { Lx, Ly, 0.0, h}; // Bottom Face +Point(104) = { 0.0, Ly, 0.0, h}; // Bottom Face + +// Base Cube +Line(101) = {101,102}; // Bottom Face +Line(102) = {102,103}; // Bottom Face +Line(103) = {103,104}; // Bottom Face +Line(104) = {104,101}; // Bottom Face + +// Base Cube +Line Loop(101) = {101:104}; + +// Base Cube +Plane Surface(101) = {101}; diff --git a/test/test_solver/test_petsc_matrix_apply_boundary.cc b/test/test_solver/test_petsc_matrix_apply_boundary.cc index 068750a25..92e598bab 100644 --- a/test/test_solver/test_petsc_matrix_apply_boundary.cc +++ b/test/test_solver/test_petsc_matrix_apply_boundary.cc @@ -1,151 +1,140 @@ /** * @file test_petsc_matrix_profile.cc * @author Aurelia Cuba Ramos * @date Wed Jul 30 12:34:08 2014 * * @brief test the applyBoundary method 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 . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "static_communicator.hh" #include "aka_common.hh" #include "aka_csr.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_utils.hh" #include "distributed_synchronizer.hh" #include "petsc_matrix.hh" #include "fe_engine.hh" #include "dof_synchronizer.hh" #include "mesh_partition_scotch.hh" using namespace akantu; int main(int argc, char *argv[]) { initialize(argc, argv); const ElementType element_type = _triangle_3; const GhostType ghost_type = _not_ghost; 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); /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ DistributedSynchronizer * communicator = NULL; if(prank == 0) { /// creation mesh mesh.read("triangle.msh"); MeshPartitionScotch * partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); delete partition; } else { communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); } FEEngine *fem = new FEEngineTemplate(mesh, spatial_dimension, "my_fem"); DOFSynchronizer dof_synchronizer(mesh, spatial_dimension); UInt nb_global_nodes = mesh.getNbGlobalNodes(); dof_synchronizer.initGlobalDOFEquationNumbers(); // fill the matrix with UInt nb_element = mesh.getNbElement(element_type); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element_type); UInt nb_dofs_per_element = spatial_dimension * nb_nodes_per_element; SparseMatrix K(nb_global_nodes * spatial_dimension, _symmetric); K.buildProfile(mesh, dof_synchronizer, spatial_dimension); Matrix element_input(nb_dofs_per_element, nb_dofs_per_element, 1); Array K_e = Array(nb_element, nb_dofs_per_element * nb_dofs_per_element, "K_e"); Array::matrix_iterator K_e_it = K_e.begin(nb_dofs_per_element, nb_dofs_per_element); Array::matrix_iterator K_e_end = K_e.end(nb_dofs_per_element, nb_dofs_per_element); for(; K_e_it != K_e_end; ++K_e_it) *K_e_it = element_input; // assemble the test matrix fem->assembleMatrix(K_e, K, spatial_dimension, element_type, ghost_type); // create petsc matrix PETScMatrix petsc_matrix(nb_global_nodes * spatial_dimension, _symmetric); petsc_matrix.buildProfile(mesh, dof_synchronizer, spatial_dimension); // add stiffness matrix to petsc matrix petsc_matrix.add(K, 1); // create boundary array: block all dofs UInt nb_nodes = mesh.getNbNodes(); Array boundary = Array(nb_nodes, spatial_dimension, true); - //boundary(2,1) = true; - - - std::cout << dof_synchronizer.getDOFGlobalID(2*2+1) < Math::getTolerance()) + return EXIT_FAILURE; + delete communicator; finalize(); - // if (std::abs(test_passed) > Math::getTolerance()) - // return EXIT_FAILURE; - - std::cout << " Test passed!!!" << std::endl; return EXIT_SUCCESS; } diff --git a/test/test_solver/test_petsc_matrix_apply_boundary.sh b/test/test_solver/test_petsc_matrix_apply_boundary.sh new file mode 100755 index 000000000..f053c4da0 --- /dev/null +++ b/test/test_solver/test_petsc_matrix_apply_boundary.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +mpirun -np 4 ./test_petsc_matrix_apply_boundary diff --git a/test/test_solver/test_petsc_matrix_profile.cc b/test/test_solver/test_petsc_matrix_diagonal.cc similarity index 61% copy from test/test_solver/test_petsc_matrix_profile.cc copy to test/test_solver/test_petsc_matrix_diagonal.cc index 1a8f2f4e7..eef514153 100644 --- a/test/test_solver/test_petsc_matrix_profile.cc +++ b/test/test_solver/test_petsc_matrix_diagonal.cc @@ -1,137 +1,146 @@ /** - * @file test_petsc_matrix_profile.cc - * @author Aurelia Cuba Ramos - * @date Wed Jul 30 12:34:08 2014 + * @file test_petsc_matrix_diagonal.cc + * @author Aurelia Isabel Cuba Ramos + * @date Wed Apr 22 09:41:14 2015 * - * @brief test the profile generation of the PETScMatrix class + * @brief test the connectivity is correctly represented in the PETScMatrix * * @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 . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "static_communicator.hh" #include "aka_common.hh" #include "aka_csr.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_utils.hh" #include "distributed_synchronizer.hh" #include "petsc_matrix.hh" #include "fe_engine.hh" #include "dof_synchronizer.hh" - +#include "dumper_paraview.hh" #include "mesh_partition_scotch.hh" using namespace akantu; int main(int argc, char *argv[]) { initialize(argc, argv); const ElementType element_type = _triangle_3; const GhostType ghost_type = _not_ghost; 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); /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ DistributedSynchronizer * communicator = NULL; if(prank == 0) { /// creation mesh mesh.read("triangle.msh"); MeshPartitionScotch * partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); delete partition; } else { communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); } + DumperParaview mesh_dumper("mesh_dumper"); + mesh_dumper.registerMesh(mesh, spatial_dimension, _not_ghost); + mesh_dumper.dump(); + /// initialize the FEEngine and the dof_synchronizer FEEngine *fem = new FEEngineTemplate(mesh, spatial_dimension, "my_fem"); DOFSynchronizer dof_synchronizer(mesh, spatial_dimension); UInt nb_global_nodes = mesh.getNbGlobalNodes(); dof_synchronizer.initGlobalDOFEquationNumbers(); - // fill the matrix with + // construct an Akantu sparse matrix, build the profile and fill the matrix for the given mesh UInt nb_element = mesh.getNbElement(element_type); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element_type); UInt nb_dofs_per_element = spatial_dimension * nb_nodes_per_element; - SparseMatrix K(nb_global_nodes * spatial_dimension, _unsymmetric); - K.buildProfile(mesh, dof_synchronizer, spatial_dimension); - Matrix element_input(nb_dofs_per_element, nb_dofs_per_element, 1); + SparseMatrix K_akantu(nb_global_nodes * spatial_dimension, _unsymmetric); + K_akantu.buildProfile(mesh, dof_synchronizer, spatial_dimension); + /// use as elemental matrices a matrix with values equal to 1 every where + Matrix element_input(nb_dofs_per_element, nb_dofs_per_element, 1.); Array K_e = Array(nb_element, nb_dofs_per_element * nb_dofs_per_element, "K_e"); Array::matrix_iterator K_e_it = K_e.begin(nb_dofs_per_element, nb_dofs_per_element); Array::matrix_iterator K_e_end = K_e.end(nb_dofs_per_element, nb_dofs_per_element); for(; K_e_it != K_e_end; ++K_e_it) *K_e_it = element_input; // assemble the test matrix - fem->assembleMatrix(K_e, K, spatial_dimension, element_type, ghost_type); - + fem->assembleMatrix(K_e, K_akantu, spatial_dimension, element_type, ghost_type); + + /// construct a PETSc matrix + PETScMatrix K_petsc(nb_global_nodes * spatial_dimension, _unsymmetric); + /// build the profile of the PETSc matrix for the mesh of this example + K_petsc.buildProfile(mesh, dof_synchronizer, spatial_dimension); + /// add an Akantu sparse matrix to a PETSc sparse matrix + K_petsc.add(K_akantu, 1); + + /// check to how many elements each node is connected CSR node_to_elem; MeshUtils::buildNode2Elements(mesh, node_to_elem, spatial_dimension); - for (UInt i = 0; i < mesh.getNbNodes(); ++i) { - std::cout << node_to_elem.getNbCols(i) << std::endl; + /// test the diagonal of the PETSc matrix: the diagonal entries + /// of the PETSc matrix correspond to the number of elements + /// connected to the node of the dof. Note: for an Akantu matrix this is only true for the serial case + Real error = 0.; + /// loop over all diagonal values of the matrix + for (UInt i = 0; i < mesh.getNbNodes(); ++i) { + for (UInt j = 0; j < spatial_dimension; ++j) { + UInt dof = i * spatial_dimension + j; + /// for PETSc matrix only DOFs on the processor and be accessed + if (dof_synchronizer.isLocalOrMasterDOF(dof)) { + UInt global_dof = dof_synchronizer.getDOFGlobalID(dof); + std::cout << "Number of elements connected: " << node_to_elem.getNbCols(i) << std::endl; + std::cout << "K_petsc(" << global_dof << "," << global_dof << ")=" << K_petsc(dof,dof) << std::endl; + error += std::abs(K_petsc(dof, dof) - node_to_elem.getNbCols(i)); + } + } } - PETScMatrix petsc_matrix(nb_global_nodes * spatial_dimension, _unsymmetric); - - petsc_matrix.buildProfile(mesh, dof_synchronizer, spatial_dimension); - - - petsc_matrix.add(K, 1); - - // for (UInt i = 0; i < mesh.getNbNodes(); ++i) { - // for (UInt j = 0; j < spatial_dimension; ++j) { - // UInt dof = i * spatial_dimension + j; - // if (dof_synchronizer.isLocalOrMasterDOF(dof)) { - // UInt global_dof = dof_synchronizer.getDOFGlobalID(dof); - // std::cout << "K(" << global_dof << "," << global_dof << ")=" << petsc_matrix(dof,dof) << std::endl; - // std::cout << node_to_elem.getNbCols(i) << std::endl; - // } - // } - // } - - - - - petsc_matrix.saveMatrix("profile.dat"); + if(error > Math::getTolerance() ) { + std::cout << "error in the stiffness matrix!!!" << std::cout; + return EXIT_FAILURE; + } delete communicator; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_solver/test_petsc_matrix_diagonal.sh b/test/test_solver/test_petsc_matrix_diagonal.sh new file mode 100755 index 000000000..5f592e4b2 --- /dev/null +++ b/test/test_solver/test_petsc_matrix_diagonal.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +mpirun -np 4 ./test_petsc_matrix_diagonal diff --git a/test/test_solver/test_petsc_matrix_profile.cc b/test/test_solver/test_petsc_matrix_profile.cc index 1a8f2f4e7..bff41dcfd 100644 --- a/test/test_solver/test_petsc_matrix_profile.cc +++ b/test/test_solver/test_petsc_matrix_profile.cc @@ -1,137 +1,131 @@ /** * @file test_petsc_matrix_profile.cc * @author Aurelia Cuba Ramos * @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 . * */ /* -------------------------------------------------------------------------- */ #include +#include /* -------------------------------------------------------------------------- */ #include "static_communicator.hh" #include "aka_common.hh" #include "aka_csr.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_utils.hh" #include "distributed_synchronizer.hh" #include "petsc_matrix.hh" #include "fe_engine.hh" #include "dof_synchronizer.hh" - +/// #include "dumper_paraview.hh" #include "mesh_partition_scotch.hh" using namespace akantu; int main(int argc, char *argv[]) { initialize(argc, argv); const ElementType element_type = _triangle_3; const GhostType ghost_type = _not_ghost; 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); /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ DistributedSynchronizer * communicator = NULL; if(prank == 0) { /// creation mesh - mesh.read("triangle.msh"); + mesh.read("square.msh"); MeshPartitionScotch * partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); delete partition; } else { communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); } - + // dump mesh in paraview + // DumperParaview mesh_dumper("mesh_dumper"); + // mesh_dumper.registerMesh(mesh, spatial_dimension, _not_ghost); + // mesh_dumper.dump(); + + /// initialize the FEEngine and the dof_synchronizer FEEngine *fem = new FEEngineTemplate(mesh, spatial_dimension, "my_fem"); DOFSynchronizer dof_synchronizer(mesh, spatial_dimension); UInt nb_global_nodes = mesh.getNbGlobalNodes(); dof_synchronizer.initGlobalDOFEquationNumbers(); - // fill the matrix with + // construct an Akantu sparse matrix, build the profile and fill the matrix for the given mesh UInt nb_element = mesh.getNbElement(element_type); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element_type); UInt nb_dofs_per_element = spatial_dimension * nb_nodes_per_element; - SparseMatrix K(nb_global_nodes * spatial_dimension, _unsymmetric); - K.buildProfile(mesh, dof_synchronizer, spatial_dimension); - Matrix element_input(nb_dofs_per_element, nb_dofs_per_element, 1); + SparseMatrix K_akantu(nb_global_nodes * spatial_dimension, _unsymmetric); + K_akantu.buildProfile(mesh, dof_synchronizer, spatial_dimension); + /// use as elemental matrices a matrix with values equal to 1 every where + Matrix element_input(nb_dofs_per_element, nb_dofs_per_element, 1.); Array K_e = Array(nb_element, nb_dofs_per_element * nb_dofs_per_element, "K_e"); Array::matrix_iterator K_e_it = K_e.begin(nb_dofs_per_element, nb_dofs_per_element); Array::matrix_iterator K_e_end = K_e.end(nb_dofs_per_element, nb_dofs_per_element); for(; K_e_it != K_e_end; ++K_e_it) *K_e_it = element_input; // assemble the test matrix - fem->assembleMatrix(K_e, K, spatial_dimension, element_type, ghost_type); - - CSR node_to_elem; - - MeshUtils::buildNode2Elements(mesh, node_to_elem, spatial_dimension); - - for (UInt i = 0; i < mesh.getNbNodes(); ++i) { - std::cout << node_to_elem.getNbCols(i) << std::endl; - } - - PETScMatrix petsc_matrix(nb_global_nodes * spatial_dimension, _unsymmetric); - - petsc_matrix.buildProfile(mesh, dof_synchronizer, spatial_dimension); - - - petsc_matrix.add(K, 1); + fem->assembleMatrix(K_e, K_akantu, spatial_dimension, element_type, ghost_type); - // for (UInt i = 0; i < mesh.getNbNodes(); ++i) { - // for (UInt j = 0; j < spatial_dimension; ++j) { - // UInt dof = i * spatial_dimension + j; - // if (dof_synchronizer.isLocalOrMasterDOF(dof)) { - // UInt global_dof = dof_synchronizer.getDOFGlobalID(dof); - // std::cout << "K(" << global_dof << "," << global_dof << ")=" << petsc_matrix(dof,dof) << std::endl; - // std::cout << node_to_elem.getNbCols(i) << std::endl; - // } - // } - // } - - - - - petsc_matrix.saveMatrix("profile.dat"); + /// construct a PETSc matrix + PETScMatrix K_petsc(nb_global_nodes * spatial_dimension, _unsymmetric); + /// build the profile of the PETSc matrix for the mesh of this example + K_petsc.buildProfile(mesh, dof_synchronizer, spatial_dimension); + /// add an Akantu sparse matrix to a PETSc sparse matrix + K_petsc.add(K_akantu, 1); + + /// save the profile + K_petsc.saveMatrix("profile.txt"); + + /// print the matrix to screen + std::ifstream profile; + profile.open("profile.txt"); + std::string current_line; + while(getline(profile, current_line)) + std::cout << current_line << std::endl; + profile.close(); delete communicator; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_solver/test_petsc_matrix_profile.verified b/test/test_solver/test_petsc_matrix_profile.verified new file mode 100644 index 000000000..4213dba00 --- /dev/null +++ b/test/test_solver/test_petsc_matrix_profile.verified @@ -0,0 +1,12 @@ +Matrix Object: 1 MPI processes + type: seqaij +row 0: (0, 2) (1, 2) (2, 1) (3, 1) (6, 1) (7, 1) (8, 2) (9, 2) +row 1: (0, 2) (1, 2) (2, 1) (3, 1) (6, 1) (7, 1) (8, 2) (9, 2) +row 2: (0, 1) (1, 1) (2, 2) (3, 2) (4, 1) (5, 1) (8, 2) (9, 2) +row 3: (0, 1) (1, 1) (2, 2) (3, 2) (4, 1) (5, 1) (8, 2) (9, 2) +row 4: (2, 1) (3, 1) (4, 2) (5, 2) (6, 1) (7, 1) (8, 2) (9, 2) +row 5: (2, 1) (3, 1) (4, 2) (5, 2) (6, 1) (7, 1) (8, 2) (9, 2) +row 6: (0, 1) (1, 1) (4, 1) (5, 1) (6, 2) (7, 2) (8, 2) (9, 2) +row 7: (0, 1) (1, 1) (4, 1) (5, 1) (6, 2) (7, 2) (8, 2) (9, 2) +row 8: (0, 2) (1, 2) (2, 2) (3, 2) (4, 2) (5, 2) (6, 2) (7, 2) (8, 4) (9, 4) +row 9: (0, 2) (1, 2) (2, 2) (3, 2) (4, 2) (5, 2) (6, 2) (7, 2) (8, 4) (9, 4) diff --git a/test/test_solver/test_petsc_matrix_profile.cc b/test/test_solver/test_petsc_matrix_profile_parallel.cc similarity index 71% copy from test/test_solver/test_petsc_matrix_profile.cc copy to test/test_solver/test_petsc_matrix_profile_parallel.cc index 1a8f2f4e7..2ec40523c 100644 --- a/test/test_solver/test_petsc_matrix_profile.cc +++ b/test/test_solver/test_petsc_matrix_profile_parallel.cc @@ -1,137 +1,130 @@ /** * @file test_petsc_matrix_profile.cc * @author Aurelia Cuba Ramos * @date Wed Jul 30 12:34:08 2014 * - * @brief test the profile generation of the PETScMatrix class + * @brief test the profile generation of the PETScMatrix class in parallel * * @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 . * */ /* -------------------------------------------------------------------------- */ #include +#include /* -------------------------------------------------------------------------- */ #include "static_communicator.hh" #include "aka_common.hh" #include "aka_csr.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_utils.hh" #include "distributed_synchronizer.hh" #include "petsc_matrix.hh" #include "fe_engine.hh" #include "dof_synchronizer.hh" - +/// #include "dumper_paraview.hh" #include "mesh_partition_scotch.hh" using namespace akantu; int main(int argc, char *argv[]) { initialize(argc, argv); const ElementType element_type = _triangle_3; const GhostType ghost_type = _not_ghost; 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); /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ DistributedSynchronizer * communicator = NULL; if(prank == 0) { /// creation mesh - mesh.read("triangle.msh"); + mesh.read("square.msh"); MeshPartitionScotch * partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); delete partition; } else { communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); } - + // dump mesh in paraview + // DumperParaview mesh_dumper("mesh_dumper"); + // mesh_dumper.registerMesh(mesh, spatial_dimension, _not_ghost); + // mesh_dumper.dump(); + + /// initialize the FEEngine and the dof_synchronizer FEEngine *fem = new FEEngineTemplate(mesh, spatial_dimension, "my_fem"); DOFSynchronizer dof_synchronizer(mesh, spatial_dimension); UInt nb_global_nodes = mesh.getNbGlobalNodes(); dof_synchronizer.initGlobalDOFEquationNumbers(); - // fill the matrix with + // construct an Akantu sparse matrix, build the profile and fill the matrix for the given mesh UInt nb_element = mesh.getNbElement(element_type); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element_type); UInt nb_dofs_per_element = spatial_dimension * nb_nodes_per_element; - SparseMatrix K(nb_global_nodes * spatial_dimension, _unsymmetric); - K.buildProfile(mesh, dof_synchronizer, spatial_dimension); - Matrix element_input(nb_dofs_per_element, nb_dofs_per_element, 1); + SparseMatrix K_akantu(nb_global_nodes * spatial_dimension, _unsymmetric); + K_akantu.buildProfile(mesh, dof_synchronizer, spatial_dimension); + /// use as elemental matrices a matrix with values equal to 1 every where + Matrix element_input(nb_dofs_per_element, nb_dofs_per_element, 1.); Array K_e = Array(nb_element, nb_dofs_per_element * nb_dofs_per_element, "K_e"); Array::matrix_iterator K_e_it = K_e.begin(nb_dofs_per_element, nb_dofs_per_element); Array::matrix_iterator K_e_end = K_e.end(nb_dofs_per_element, nb_dofs_per_element); for(; K_e_it != K_e_end; ++K_e_it) *K_e_it = element_input; // assemble the test matrix - fem->assembleMatrix(K_e, K, spatial_dimension, element_type, ghost_type); - - CSR node_to_elem; - - MeshUtils::buildNode2Elements(mesh, node_to_elem, spatial_dimension); - - for (UInt i = 0; i < mesh.getNbNodes(); ++i) { - std::cout << node_to_elem.getNbCols(i) << std::endl; - } - - PETScMatrix petsc_matrix(nb_global_nodes * spatial_dimension, _unsymmetric); - - petsc_matrix.buildProfile(mesh, dof_synchronizer, spatial_dimension); - - - petsc_matrix.add(K, 1); + fem->assembleMatrix(K_e, K_akantu, spatial_dimension, element_type, ghost_type); - // for (UInt i = 0; i < mesh.getNbNodes(); ++i) { - // for (UInt j = 0; j < spatial_dimension; ++j) { - // UInt dof = i * spatial_dimension + j; - // if (dof_synchronizer.isLocalOrMasterDOF(dof)) { - // UInt global_dof = dof_synchronizer.getDOFGlobalID(dof); - // std::cout << "K(" << global_dof << "," << global_dof << ")=" << petsc_matrix(dof,dof) << std::endl; - // std::cout << node_to_elem.getNbCols(i) << std::endl; - // } - // } - // } - - - - - petsc_matrix.saveMatrix("profile.dat"); + /// construct a PETSc matrix + PETScMatrix K_petsc(nb_global_nodes * spatial_dimension, _unsymmetric); + /// build the profile of the PETSc matrix for the mesh of this example + K_petsc.buildProfile(mesh, dof_synchronizer, spatial_dimension); + /// add an Akantu sparse matrix to a PETSc sparse matrix + K_petsc.add(K_akantu, 1); + + /// save the profile + K_petsc.saveMatrix("profile_parallel.txt"); + /// print the matrix to screen + std::ifstream profile; + profile.open("profile_parallel.txt"); + std::string current_line; + while(getline(profile, current_line)) + std::cout << current_line << std::endl; + profile.close(); delete communicator; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_solver/test_petsc_matrix_profile_parallel.sh b/test/test_solver/test_petsc_matrix_profile_parallel.sh new file mode 100755 index 000000000..2dd775687 --- /dev/null +++ b/test/test_solver/test_petsc_matrix_profile_parallel.sh @@ -0,0 +1,3 @@ +#!/bin/bash + +mpirun -np 2 ./test_petsc_matrix_profile_parallel diff --git a/test/test_solver/test_petsc_matrix_profile_parallel.verified b/test/test_solver/test_petsc_matrix_profile_parallel.verified new file mode 100644 index 000000000..662687ada --- /dev/null +++ b/test/test_solver/test_petsc_matrix_profile_parallel.verified @@ -0,0 +1,12 @@ +Matrix Object: 1 MPI processes + type: mpiaij +row 0: (0, 2) (1, 2) (2, 1) (3, 1) (4, 1) (5, 1) (6, 2) (7, 2) +row 1: (0, 2) (1, 2) (2, 1) (3, 1) (4, 1) (5, 1) (6, 2) (7, 2) +row 2: (0, 1) (1, 1) (2, 2) (3, 2) (6, 2) (7, 2) (8, 1) (9, 1) +row 3: (0, 1) (1, 1) (2, 2) (3, 2) (6, 2) (7, 2) (8, 1) (9, 1) +row 4: (0, 1) (1, 1) (4, 2) (5, 2) (6, 2) (7, 2) (8, 1) (9, 1) +row 5: (0, 1) (1, 1) (4, 2) (5, 2) (6, 2) (7, 2) (8, 1) (9, 1) +row 6: (0, 2) (1, 2) (2, 2) (3, 2) (4, 2) (5, 2) (6, 4) (7, 4) (8, 2) (9, 2) +row 7: (0, 2) (1, 2) (2, 2) (3, 2) (4, 2) (5, 2) (6, 4) (7, 4) (8, 2) (9, 2) +row 8: (2, 1) (3, 1) (4, 1) (5, 1) (6, 2) (7, 2) (8, 2) (9, 2) +row 9: (2, 1) (3, 1) (4, 1) (5, 1) (6, 2) (7, 2) (8, 2) (9, 2) diff --git a/test/test_solver/test_solver_petsc.cc b/test/test_solver/test_solver_petsc.cc index 41aaea640..2eb172f50 100644 --- a/test/test_solver/test_solver_petsc.cc +++ b/test/test_solver/test_solver_petsc.cc @@ -1,154 +1,161 @@ /** * @file test_solver_petsc.cc * @author Aurelia Cuba Ramos * @date Tue Dec 2 17:17:34 2014 * * @brief test the PETSc solver interface * * @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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "static_communicator.hh" #include "aka_common.hh" #include "aka_csr.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_utils.hh" #include "distributed_synchronizer.hh" #include "petsc_matrix.hh" #include "solver_petsc.hh" #include "fe_engine.hh" #include "dof_synchronizer.hh" #include "mesh_partition_scotch.hh" using namespace akantu; int main(int argc, char *argv[]) { initialize(argc, argv); const ElementType element_type = _segment_2; const GhostType ghost_type = _not_ghost; UInt spatial_dimension = 1; StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); /// read the mesh and partition it Mesh mesh(spatial_dimension); /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ DistributedSynchronizer * communicator = NULL; if(prank == 0) { /// creation mesh mesh.read("1D_bar.msh"); MeshPartitionScotch * partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); delete partition; } else { communicator = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); } FEEngine *fem = new FEEngineTemplate(mesh, spatial_dimension, "my_fem"); DOFSynchronizer dof_synchronizer(mesh, spatial_dimension); UInt nb_global_nodes = mesh.getNbGlobalNodes(); dof_synchronizer.initGlobalDOFEquationNumbers(); // fill the matrix with UInt nb_element = mesh.getNbElement(element_type); - std::cout << mesh.getNbElement(element_type, _ghost) < element_input(nb_dofs_per_element, nb_dofs_per_element, 0); for (UInt i = 0; i < nb_dofs_per_element; ++i) { for (UInt j = 0; j < nb_dofs_per_element; ++j) { element_input(i, j) = ((i == j) ? 1 : -1); } } Array K_e = Array(nb_element, nb_dofs_per_element * nb_dofs_per_element, "K_e"); Array::matrix_iterator K_e_it = K_e.begin(nb_dofs_per_element, nb_dofs_per_element); Array::matrix_iterator K_e_end = K_e.end(nb_dofs_per_element, nb_dofs_per_element); for(; K_e_it != K_e_end; ++K_e_it) *K_e_it = element_input; // assemble the test matrix fem->assembleMatrix(K_e, K, spatial_dimension, element_type, ghost_type); // apply boundary: block first node const Array & position = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); Array boundary = Array(nb_nodes, spatial_dimension, false); for (UInt i = 0; i < nb_nodes; ++i) { if (std::abs(position(i, 0)) < Math::getTolerance() ) boundary(i, 0) = true; } K.applyBoundary(boundary); /// create the PETSc matrix for the solve step PETScMatrix petsc_matrix(nb_global_nodes * spatial_dimension, _symmetric); petsc_matrix.buildProfile(mesh, dof_synchronizer, spatial_dimension); /// copy the stiffness matrix into the petsc matrix petsc_matrix.add(K, 1); // initialize internal forces: they are zero because imposed displacement is zero Array internal_forces(nb_nodes, spatial_dimension, 0.); // compute residual: apply nodal force on last node Array residual(nb_nodes, spatial_dimension, 0.); - residual -= internal_forces; + for (UInt i = 0; i < nb_nodes; ++i) { if (std::abs(position(i, 0) - 10) < Math::getTolerance() ) residual(i, 0) += 2; } - /// solve the matrix before the solve step - petsc_matrix.saveMatrix("1D_bar.mtx"); + residual -= internal_forces; /// initialize solver and solution Array solution(nb_nodes, spatial_dimension, 0.); SolverPETSc solver(petsc_matrix); solver.initialize(); + solver.setOperators(); solver.setRHS(residual); solver.solve(solution); + /// verify solution + Math::setTolerance(1e-11); + for (UInt i = 0; i < nb_nodes; ++i) { + if (!dof_synchronizer.isPureGhostDOF(i) && !Math::are_float_equal(2 * position(i, 0), solution(i, 0))) { + std::cout << "The solution is not correct!!!!" << std::endl; + return EXIT_FAILURE; + } + } + delete communicator; finalize(); return EXIT_SUCCESS; } diff --git a/test/test_solver/test_solver_petsc.sh b/test/test_solver/test_solver_petsc.sh new file mode 100755 index 000000000..edcc7916c --- /dev/null +++ b/test/test_solver/test_solver_petsc.sh @@ -0,0 +1,5 @@ +#! /bin/bash + +# choose solver Mumps through the PETSc interface +./test_solver_petsc -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps -mat_mumps_icntl_7 2 +mpirun -np 4 ./test_solver_petsc -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps -mat_mumps_icntl_7 2