diff --git a/packages/core.cmake b/packages/core.cmake
index 9b1865be2..2ecc7a25a 100644
--- a/packages/core.cmake
+++ b/packages/core.cmake
@@ -1,234 +1,221 @@
#===============================================================================
# @file core.cmake
#
# @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date Mon Nov 21 18:19:15 2011
#
# @brief package description for core
#
# @section LICENSE
#
# Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# Akantu is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Akantu. If not, see <http://www.gnu.org/licenses/>.
#
#===============================================================================
set(AKANTU_CORE ON CACHE INTERNAL "core package for Akantu" FORCE)
set(AKANTU_CORE_FILES
# source files
common/aka_common.cc
common/aka_error.cc
common/aka_extern.cc
common/aka_static_memory.cc
common/aka_memory.cc
common/aka_vector.cc
common/aka_math.cc
fem/shape_lagrange.cc
fem/shape_linked.cc
- #fem/integrator_gauss.cc
fem/mesh.cc
fem/fem.cc
io/dumper/dumpable.hh
io/mesh_io.cc
io/model_io.cc
io/mesh_io/mesh_io_msh.cc
io/mesh_io/mesh_io_diana.cc
model/model.cc
model/solid_mechanics/solid_mechanics_model.cc
model/solid_mechanics/solid_mechanics_model_mass.cc
model/solid_mechanics/solid_mechanics_model_boundary.cc
model/solid_mechanics/solid_mechanics_model_material.cc
model/solid_mechanics/material.cc
model/solid_mechanics/material_parameters.cc
model/solid_mechanics/materials/material_elastic.cc
mesh_utils/mesh_pbc.cc
mesh_utils/mesh_partition.cc
mesh_utils/mesh_utils.cc
solver/sparse_matrix.cc
solver/solver.cc
synchronizer/synchronizer_registry.cc
synchronizer/synchronizer.cc
synchronizer/distributed_synchronizer.cc
synchronizer/pbc_synchronizer.cc
synchronizer/data_accessor.cc
synchronizer/static_communicator.cc
synchronizer/dof_synchronizer.cc
#header files
io/mesh_io.hh
io/model_io.hh
io/mesh_io/mesh_io_msh.hh
io/mesh_io/mesh_io_diana.hh
mesh_utils/mesh_utils.hh
mesh_utils/mesh_partition.hh
mesh_utils/mesh_partition/mesh_partition_scotch.hh
solver/sparse_matrix.hh
solver/solver.hh
synchronizer/synchronizer.hh
synchronizer/synchronizer_registry.hh
synchronizer/static_communicator_dummy.hh
synchronizer/static_communicator_inline_impl.hh
synchronizer/distributed_synchronizer.hh
synchronizer/pbc_synchronizer.hh
synchronizer/static_communicator.hh
synchronizer/dof_synchronizer.hh
synchronizer/real_static_communicator.hh
synchronizer/data_accessor.hh
synchronizer/communication_buffer.hh
common/aka_fwd.hh
common/aka_grid.hh
common/aka_grid_tmpl.hh
common/aka_types.hh
common/aka_static_memory.hh
common/aka_static_memory_tmpl.hh
common/aka_memory.hh
common/aka_math.hh
common/aka_math_tmpl.hh
common/aka_csr.hh
common/aka_error.hh
common/aka_common.hh
common/aka_vector.hh
common/aka_vector_tmpl.hh
common/aka_circular_vector.hh
common/aka_event_handler.hh
common/aka_random_generator.hh
common/aka_bounding_box.hh
common/aka_ci_string.hh
common/aka_plane.hh
common/aka_polytope.hh
common/aka_sphere.hh
common/aka_timer.hh
common/aka_tree.hh
common/aka_typelist.hh
common/aka_visitor.hh
+ common/aka_grid_dynamic.hh
+ common/aka_safe_enum.hh
fem/mesh.hh
fem/fem.hh
fem/by_element_type.hh
fem/shape_functions.hh
fem/shape_lagrange.hh
fem/fem_template.hh
fem/fem_template_tmpl.hh
fem/integrator_gauss.hh
fem/integrator.hh
fem/element_class.hh
fem/shape_linked.hh
+ fem/geometrical_data_tmpl.hh
model/model.hh
model/parser.hh
model/parser_tmpl.hh
model/structural_mechanics/structural_mechanics_model.hh
model/integration_scheme/integration_scheme_2nd_order.hh
model/integration_scheme/generalized_trapezoidal.hh
model/integration_scheme/newmark-beta.hh
model/integration_scheme/integration_scheme_1st_order.hh
model/solid_mechanics/solid_mechanics_model.hh
model/solid_mechanics/solid_mechanics_model_tmpl.hh
model/solid_mechanics/material.hh
model/solid_mechanics/material_parameters.hh
model/solid_mechanics/material_parameters_tmpl.hh
model/solid_mechanics/materials/material_elastic.hh
#inline implementation files
mesh_utils/mesh_utils_inline_impl.cc
solver/sparse_matrix_inline_impl.cc
synchronizer/dof_synchronizer_inline_impl.cc
synchronizer/communication_buffer_inline_impl.cc
common/aka_common_inline_impl.cc
common/aka_memory_inline_impl.cc
common/aka_static_memory_inline_impl.cc
common/aka_circular_vector_inline_impl.cc
fem/integrator_gauss_inline_impl.cc
fem/fem_template_inline_impl.cc
fem/element_classes/element_class_triangle_3_inline_impl.cc
fem/element_classes/element_class_segment_2_inline_impl.cc
fem/element_classes/element_class_quadrangle_4_inline_impl.cc
fem/element_classes/element_class_quadrangle_8_inline_impl.cc
fem/element_classes/element_class_hexahedron_8_inline_impl.cc
fem/element_classes/element_class_triangle_6_inline_impl.cc
fem/element_classes/element_class_tetrahedron_10_inline_impl.cc
fem/element_classes/element_class_segment_3_inline_impl.cc
fem/element_classes/element_class_tetrahedron_4_inline_impl.cc
fem/shape_functions_inline_impl.cc
fem/mesh_inline_impl.cc
fem/by_element_type_tmpl.hh
fem/fem_inline_impl.cc
fem/shape_linked_inline_impl.cc
fem/shape_lagrange_inline_impl.cc
model/model_inline_impl.cc
model/integration_scheme/generalized_trapezoidal_inline_impl.cc
model/integration_scheme/newmark-beta_inline_impl.cc
model/solid_mechanics/solid_mechanics_model_inline_impl.cc
model/solid_mechanics/materials/material_elastic_inline_impl.cc
model/solid_mechanics/material_inline_impl.cc
model/parser_inline_impl.cc
fem/geometrical_element.cc
fem/element_class_tmpl.hh
fem/element_class.cc
fem/integration_element.cc
fem/interpolation_element.cc
model/solid_mechanics/materials/material_elastic_orthotropic.cc
model/solid_mechanics/materials/material_elastic_orthotropic.hh
model/solid_mechanics/materials/material_elastic_orthotropic_inline_impl.cc
)
-#include(CheckCXXCompilerFlag)
-#check_cxx_compiler_flag (-std=c++0x HAVE_NEW_STD)
-#if (HAVE_NEW_STD)
-# list(APPEND AKANTU_CORE_FILES
-# common/aka_point.hh
-# common/aka_bounding_box.hh
-# common/aka_bounding_box.cc
-# common/aka_geometry.hh
-# common/aka_geometry.cc
-# )
-# add_definitions(-std=c++0x)
-#else()
-# message(WARNING "*** WARNING *** Compiler does not support c++11 set of requirements.")
-#endif()
-
set(AKANTU_CORE_DEB_DEPEND
libboost-dev
)
set(AKANTU_CORE_TESTS
test_solid_mechanics_model_square
test_vector
test_vector_iterator
test_matrix
test_csr
test_grid
test_static_memory
test_mesh_io_msh
test_facet_extraction_triangle_3
test_facet_extraction_tetrahedron_4
test_pbc_tweak
test_purify_mesh
test_local_material
test_interpolate_stress
test_weight
test_solid_mechanics_model_circle_2
test_solid_mechanics_model_bar_traction2d
test_solid_mechanics_model_bar_traction2d_structured
test_solid_mechanics_model_bar_traction2d_structured_pbc
test_solid_mechanics_model_cube3d
test_solid_mechanics_model_cube3d_tetra10
test_solid_mechanics_model_cube3d_pbc
test_surface_extraction_triangle_3
test_surface_extraction_tetrahedron_4
test_material_damage_non_local
)
\ No newline at end of file
diff --git a/src/common/aka_grid_dynamic.hh b/src/common/aka_grid_dynamic.hh
index 0b85b5768..808b07dbf 100644
--- a/src/common/aka_grid_dynamic.hh
+++ b/src/common/aka_grid_dynamic.hh
@@ -1,451 +1,452 @@
/**
* @file aka_grid_dynamic.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Feb 15 16:55:32 2013
*
* @brief Grid that is auto balanced
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_types.hh"
#include <iostream>
/* -------------------------------------------------------------------------- */
#include <map>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_GRID_DYNAMIC_HH__
#define __AKANTU_AKA_GRID_DYNAMIC_HH__
__BEGIN_AKANTU__
class Mesh;
template<typename T>
class SpatialGrid {
public:
SpatialGrid(UInt dimension) : dimension(dimension),
spacing(dimension),
center(dimension),
lower(dimension),
upper(dimension),
empty_cell() {}
SpatialGrid(UInt dimension,
const types::Vector<Real> & spacing,
const types::Vector<Real> & center) : dimension(dimension),
spacing(spacing),
center(center),
lower(dimension),
upper(dimension),
empty_cell() {
for (UInt i = 0; i < dimension; ++i) {
lower(i) = std::numeric_limits<Real>::max();
upper(i) = - std::numeric_limits<Real>::max();
}
}
virtual ~SpatialGrid() {};
class neighbor_cells_iterator;
class CellID {
public:
CellID() : ids() {}
CellID(UInt dimention) : ids(dimention) {}
void setID(UInt dir, Int id) { ids(dir) = id; }
Int getID(UInt dir) const { return ids(dir); }
bool operator<(const CellID & id) const {
return std::lexicographical_compare(ids.storage(), ids.storage() + ids.size(),
id.ids.storage(), id.ids.storage() + id.ids.size());
}
bool operator==(const CellID & id) const {
return std::equal(ids.storage(), ids.storage() + ids.size(),
id.ids.storage());
}
bool operator!=(const CellID & id) const {
return !(operator==(id));
}
private:
friend class neighbor_cells_iterator;
types::Vector<Int> ids;
};
/* -------------------------------------------------------------------------- */
class Cell {
public:
typedef typename std::vector<T>::iterator iterator;
typedef typename std::vector<T>::const_iterator const_iterator;
Cell() : id(), data() { }
Cell(const CellID &cell_id) : id(cell_id), data() { }
bool operator==(const Cell & cell) const { return id == cell.id; }
bool operator!=(const Cell & cell) const { return id != cell.id; }
Cell & add(const T & d) { data.push_back(d); return *this; }
iterator begin() { return data.begin(); }
const_iterator begin() const { return data.begin(); }
iterator end() { return data.end(); }
const_iterator end() const { return data.end(); }
#if not defined(AKANTU_NDEBUG)
Cell & add(const T & d, const types::Vector<Real> & pos) {
data.push_back(d); positions.push_back(pos); return *this;
}
typedef typename std::vector< types::Vector<Real> >::const_iterator position_iterator;
position_iterator begin_pos() const { return positions.begin(); }
position_iterator end_pos() const { return positions.end(); }
#endif
private:
CellID id;
std::vector<T> data;
#if not defined(AKANTU_NDEBUG)
std::vector< types::Vector<Real> > positions;
#endif
};
private:
typedef std::map<CellID, Cell> cells_container;
public:
const Cell & getCell(const CellID & cell_id) const {
typename cells_container::const_iterator it = cells.find(cell_id);
if(it != cells.end()) return it->second;
else return empty_cell;
}
typename Cell::iterator beginCell(const CellID & cell_id) {
typename cells_container::iterator it = cells.find(cell_id);
if(it != cells.end()) return it->second.begin();
else return empty_cell.begin();
}
typename Cell::iterator endCell(const CellID & cell_id) {
typename cells_container::iterator it = cells.find(cell_id);
if(it != cells.end()) return it->second.end();
else return empty_cell.end();
}
typename Cell::const_iterator beginCell(const CellID & cell_id) const {
typename cells_container::const_iterator it = cells.find(cell_id);
if(it != cells.end()) return it->second.begin();
else return empty_cell.begin();
}
typename Cell::const_iterator endCell(const CellID & cell_id) const {
typename cells_container::const_iterator it = cells.find(cell_id);
if(it != cells.end()) return it->second.end();
else return empty_cell.end();
}
class neighbor_cells_iterator : private std::iterator<std::forward_iterator_tag, UInt> {
public:
neighbor_cells_iterator(const CellID & cell_id, bool end) :
cell_id(cell_id),
position(cell_id.ids.size(), end ? 1 : -1) {
this->updateIt();
if(end) this->it++;
}
neighbor_cells_iterator& operator++() {
UInt i = 0;
for (; i < position.size() && position(i) == 1; ++i);
if(i == position.size()) ++it;
else {
for (UInt j = 0; j < i; ++j) position(j) = -1;
position(i)++;
updateIt();
}
return *this;
}
neighbor_cells_iterator operator++(int) { neighbor_cells_iterator tmp(*this); operator++(); return tmp; };
bool operator==(const neighbor_cells_iterator& rhs) const { return cell_id == rhs.cell_id && it == rhs.it; };
bool operator!=(const neighbor_cells_iterator& rhs) const { return ! operator==(rhs); };
CellID operator*() const {
CellID cur_cell_id(cell_id);
cur_cell_id.ids += position;
return cur_cell_id;
};
private:
void updateIt() {
it = 0;
for (UInt i = 0; i < position.size(); ++i) it = it * 3 + (position(i) + 1);
}
private:
/// central cell id
const CellID & cell_id;
// number representing the current neighbor in base 3;
UInt it;
types::Vector<Int> position;
};
public:
Cell & insert(const T & d, const types::Vector<Real> & position) {
CellID cell_id = getCellID(position);
typename cells_container::iterator it = cells.find(cell_id);
if(it == cells.end()) {
Cell cell(cell_id);
#if defined(AKANTU_NDEBUG)
Cell & tmp = (cells[cell_id] = cell).add(d);
#else
Cell & tmp = (cells[cell_id] = cell).add(d, position);
#endif
for (UInt i = 0; i < dimension; ++i) {
- Real posl = cell_id.getID(i) * spacing(i);
+ Real posl = center(i) + cell_id.getID(i) * spacing(i);
Real posu = posl + spacing(i);
if(posl < lower(i)) lower(i) = posl;
if(posu > upper(i)) upper(i) = posu;
}
return tmp;
} else {
#if defined(AKANTU_NDEBUG)
return it->second.add(d);
#else
return it->second.add(d, position);
#endif
}
}
inline neighbor_cells_iterator beginNeighborCells(const CellID & cell_id) const {
return neighbor_cells_iterator(cell_id, false);
}
inline neighbor_cells_iterator endNeighborCells(const CellID & cell_id) const {
return neighbor_cells_iterator(cell_id, true);
}
- CellID getCellID(types::Vector<Real> position) {
+ CellID getCellID(types::Vector<Real> position) const {
CellID cell_id(dimension);
for (UInt i = 0; i < dimension; ++i) {
cell_id.setID(i, getCellID(position(i), i));
}
return cell_id;
}
void printself(std::ostream & stream, int indent = 0) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
std::streamsize prec = stream.precision();
std::ios_base::fmtflags ff = stream.flags();
stream.setf (std::ios_base::showbase);
stream.precision(5);
stream << space << "SpatialGrid<" << debug::demangle(typeid(T).name()) << "> [" << std::endl;
stream << space << " + dimension : " << this->dimension << std::endl;
stream << space << " + lower bounds : {";
for (UInt i = 0; i < lower.size(); ++i) { if(i != 0) stream << ", "; stream << lower(i); };
stream << "}" << std::endl;
stream << space << " + upper bounds : {";
for (UInt i = 0; i < upper.size(); ++i) { if(i != 0) stream << ", "; stream << upper(i); };
stream << "}" << std::endl;
stream << space << " + spacing : {";
for (UInt i = 0; i < spacing.size(); ++i) { if(i != 0) stream << ", "; stream << spacing(i); };
stream << "}" << std::endl;
stream << space << " + center : {";
for (UInt i = 0; i < center.size(); ++i) { if(i != 0) stream << ", "; stream << center(i); };
stream << "}" << std::endl;
stream << space << " + nb_cells : " << this->cells.size() << "/";
types::Vector<Real> dist(this->dimension);
dist = upper;
dist -= lower;
dist /= spacing;
UInt nb_cells = std::ceil(dist(0));
for (UInt i = 1; i < this->dimension; ++i) {
nb_cells *= std::ceil(dist(i));
}
stream << nb_cells << std::endl;
stream << space << "]" << std::endl;
stream.precision(prec);
stream.flags(ff);
}
void saveAsMesh(Mesh & mesh) const;
private:
/* -------------------------------------------------------------------------- */
inline UInt getCellID(Real position, UInt direction) const {
AKANTU_DEBUG_ASSERT(direction < center.size(), "The direction asked ("
<< direction << ") is out of range " << center.size());
Real dist_center = position - center(direction);
Int id = std::floor(dist_center / spacing(direction));
//if(dist_center < 0) id--;
return id;
}
friend class GridSynchronizer;
public:
AKANTU_GET_MACRO(LowerBounds, lower, const types::Vector<Real> &);
AKANTU_GET_MACRO(UpperBounds, upper, const types::Vector<Real> &);
AKANTU_GET_MACRO(Spacing, spacing, const types::Vector<Real> &);
protected:
UInt dimension;
cells_container cells;
types::Vector<Real> spacing;
types::Vector<Real> center;
types::Vector<Real> lower;
types::Vector<Real> upper;
Cell empty_cell;
};
/// standard output stream operator
template<typename T>
inline std::ostream & operator <<(std::ostream & stream, const SpatialGrid<T> & _this)
{
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
#include "mesh.hh"
__BEGIN_AKANTU__
template<typename T>
void SpatialGrid<T>::saveAsMesh(Mesh & mesh) const {
Vector<Real> & nodes = const_cast<Vector<Real> &>(mesh.getNodes());
ElementType type;
switch(dimension) {
case 1: type = _segment_2; break;
case 2: type = _quadrangle_4; break;
case 3: type = _hexahedron_8; break;
}
+
mesh.addConnectivityType(type);
Vector<UInt> & connectivity = const_cast<Vector<UInt> &>(mesh.getConnectivity(type));
Vector<UInt> & uint_data = *mesh.getUIntDataPointer(type, "tag_1");
typename cells_container::const_iterator it = cells.begin();
typename cells_container::const_iterator end = cells.end();
types::Vector<Real> pos(dimension);
UInt global_id = 0;
for (;it != end; ++it, ++global_id) {
UInt cur_node = nodes.getSize();
UInt cur_elem = connectivity.getSize();
const CellID & cell_id = it->first;
for (UInt i = 0; i < dimension; ++i) pos(i) = center(i) + cell_id.getID(i) * spacing(i);
nodes.push_back(pos);
for (UInt i = 0; i < dimension; ++i) pos(i) += spacing(i);
nodes.push_back(pos);
connectivity.push_back(cur_node);
switch(dimension) {
case 1:
connectivity(cur_elem, 1) = cur_node + 1;
break;
case 2:
pos(0) -= spacing(0); nodes.push_back(pos);
pos(0) += spacing(0); pos(1) -= spacing(1); nodes.push_back(pos);
connectivity(cur_elem, 1) = cur_node + 3;
connectivity(cur_elem, 2) = cur_node + 1;
connectivity(cur_elem, 3) = cur_node + 2;
break;
case 3:
pos(1) -= spacing(1); pos(2) -= spacing(2); nodes.push_back(pos);
pos(1) += spacing(1); nodes.push_back(pos);
pos(0) -= spacing(0); nodes.push_back(pos);
pos(1) -= spacing(1); pos(2) += spacing(2); nodes.push_back(pos);
pos(0) += spacing(0); nodes.push_back(pos);
pos(0) -= spacing(0); pos(1) += spacing(1); nodes.push_back(pos);
connectivity(cur_elem, 1) = cur_node + 2;
connectivity(cur_elem, 2) = cur_node + 3;
connectivity(cur_elem, 3) = cur_node + 4;
connectivity(cur_elem, 4) = cur_node + 5;
connectivity(cur_elem, 5) = cur_node + 6;
connectivity(cur_elem, 6) = cur_node + 1;
connectivity(cur_elem, 7) = cur_node + 7;
break;
}
uint_data.push_back(global_id);
}
#if not defined(AKANTU_NDEBUG)
mesh.addConnectivityType(_point_1);
Vector<UInt> & connectivity_pos = const_cast<Vector<UInt> &>(mesh.getConnectivity(_point_1));
Vector<UInt> & uint_data_pos = *mesh.getUIntDataPointer(_point_1, "tag_1");
Vector<UInt> & uint_data_pos_ghost = *mesh.getUIntDataPointer(_point_1, "tag_0");
it = cells.begin();
global_id = 0;
for (;it != end; ++it, ++global_id) {
typename Cell::position_iterator cell_it = it->second.begin_pos();
typename Cell::const_iterator cell_it_cont = it->second.begin();
typename Cell::position_iterator cell_end = it->second.end_pos();
for (;cell_it != cell_end; ++cell_it, ++cell_it_cont) {
nodes.push_back(*cell_it);
connectivity_pos.push_back(nodes.getSize()-1);
uint_data_pos.push_back(global_id);
uint_data_pos_ghost.push_back(cell_it_cont->ghost_type==_ghost);
}
}
#endif
}
__END_AKANTU__
#endif /* __AKANTU_AKA_GRID_DYNAMIC_HH__ */
diff --git a/src/common/aka_math_tmpl.hh b/src/common/aka_math_tmpl.hh
index 0a4d37d90..745f42a31 100644
--- a/src/common/aka_math_tmpl.hh
+++ b/src/common/aka_math_tmpl.hh
@@ -1,728 +1,728 @@
/**
* @file aka_math_tmpl.hh
*
* @author Leonardo Snozzi <leonardo.snozzi@epfl.ch>
* @author Peter Spijker <peter.spijker@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Alejandro Marcos Aragon <alejandro.aragon@epfl.ch>
* @author Mathilde Radiguet <mathilde.radiguet@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
*
* @date Wed Aug 04 10:58:42 2010
*
* @brief Implementation of the inline functions of the math toolkit
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
__END_AKANTU__
#include <cmath>
#include <cstring>
#include <typeinfo>
#ifdef AKANTU_USE_BLAS
# ifndef AKANTU_USE_BLAS_MKL
# include <cblas.h>
# else // AKANTU_USE_BLAS_MKL
# include <mkl_cblas.h>
# endif //AKANTU_USE_BLAS_MKL
#endif //AKANTU_USE_BLAS
#ifdef AKANTU_USE_LAPACK
extern "C" {
void dgeev_(char* jobvl, char* jobvr, int* n, double* a,
int* lda, double* wr, double* wi, double* vl, int* ldvl,
double* vr, int* ldvr, double* work, int* lwork, int* info);
// LU decomposition of a general matrix
void dgetrf_(int* m, int *n,
double* a, int* lda,
int* ipiv, int* info);
// generate inverse of a matrix given its LU decomposition
void dgetri_(int* n, double* a, int* lda,
int* ipiv, double* work, int* lwork, int* info);
}
#endif
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
inline void Math::matrix_vector(UInt m, UInt n,
const Real * A,
const Real * x,
Real * y, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// y = alpha*op(A)*x + beta*y
cblas_dgemv(CblasColMajor, CblasNoTrans,
m, n, alpha, A, n, x, 1, 0, y, 1);
#else
memset(y, 0, m*sizeof(Real));
for (UInt i = 0; i < m; ++i) {
for (UInt j = 0; j < n; ++j) {
y[i] += A[i + j*m] * x[j];
}
y[i] *= alpha;
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrixt_vector(UInt m, UInt n,
const Real * A,
const Real * x,
Real * y, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// y = alpha*op(A)*x + beta*y
cblas_dgemv(CblasColMajor, CblasTrans,
m, n, alpha, A, m, x, 1, 0, y, 1);
#else
memset(y, 0, m*sizeof(Real));
for (UInt i = 0; i < m; ++i) {
for (UInt j = 0; j < n; ++j) {
y[i] += A[i * n + j] * x[j];
}
y[i] *= alpha;
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrix_matrix(UInt m, UInt n, UInt k,
const Real * A,
const Real * B,
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
m, n, k,
alpha,
A, k,
B, n,
0,
C, n);
#else
memset(C, 0, m*n*sizeof(Real));
for (UInt j = 0; j < n; ++j) {
UInt _jb = j * k;
UInt _jc = j * m;
for (UInt i = 0; i < m; ++i) {
for (UInt l = 0; l < k; ++l) {
UInt _la = l * m;
C[i + _jc] += A[i + _la] * B[l + _jb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrixt_matrix(UInt m, UInt n, UInt k,
const Real * A,
const Real * B,
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans,
m, n, k,
alpha,
A, m,
B, n,
0,
C, n);
#else
memset(C, 0, m*n*sizeof(Real));
for (UInt j = 0; j < n; ++j) {
UInt _jc = j*m;
UInt _jb = j*k;
for (UInt i = 0; i < m; ++i) {
UInt _ia = i*k;
for (UInt l = 0; l < k; ++l) {
C[i + _jc] += A[l + _ia] * B[l + _jb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrix_matrixt(UInt m, UInt n, UInt k,
const Real * A,
const Real * B,
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
m, n, k,
alpha,
A, k,
B, k,
0,
C, n);
#else
memset(C, 0, m*n*sizeof(Real));
for (UInt j = 0; j < n; ++j) {
UInt _jc = j * m;
for (UInt i = 0; i < m; ++i) {
for (UInt l = 0; l < k; ++l) {
UInt _la = l * m;
UInt _lb = l * n;
C[i + _jc] += A[i + _la] * B[j + _lb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* -------------------------------------------------------------------------- */
inline void Math::matrixt_matrixt(UInt m, UInt n, UInt k,
const Real * A,
const Real * B,
Real * C, Real alpha) {
#ifdef AKANTU_USE_BLAS
/// C := alpha*op(A)*op(B) + beta*C
cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans,
m, n, k,
alpha,
A, m,
B, k,
0,
C, n);
#else
memset(C, 0, m * n * sizeof(Real));
for (UInt j = 0; j < n; ++j) {
UInt _jc = j*m;
for (UInt i = 0; i < m; ++i) {
UInt _ia = i*k;
for (UInt l = 0; l < k; ++l) {
UInt _lb = l * n;
C[i + _jc] += A[l + _ia] * B[j + _lb];
}
C[i + _jc] *= alpha;
}
}
#endif
}
/* -------------------------------------------------------------------------- */
inline Real Math::vectorDot(const Real * v1, const Real * v2, UInt n) {
#ifdef AKANTU_USE_BLAS
/// d := v1 . v2
Real d = cblas_ddot(n, v1, 1, v2, 1);
#else
Real d = 0;
for (UInt i = 0; i < n; ++i) {
d += v1[i] * v2[i];
}
#endif
return d;
}
/* -------------------------------------------------------------------------- */
template <bool tr_A, bool tr_B>
inline void Math::matMul(UInt m, UInt n, UInt k,
Real alpha, const Real * A, const Real * B,
__attribute__ ((unused)) Real beta, Real * C) {
if(tr_A) {
if(tr_B) matrixt_matrixt(m, n, k, A, B, C, alpha);
else matrixt_matrix(m, n, k, A, B, C, alpha);
} else {
if(tr_B) matrix_matrixt(m, n, k, A, B, C, alpha);
else matrix_matrix(m, n, k, A, B, C, alpha);
}
}
/* -------------------------------------------------------------------------- */
template <bool tr_A>
inline void Math::matVectMul(UInt m, UInt n,
Real alpha, const Real * A, const Real * x,
__attribute__ ((unused)) Real beta, Real * y) {
if(tr_A) {
matrixt_vector(m, n, A, x, y, alpha);
} else {
matrix_vector(m, n, A, x, y, alpha);
}
}
/* -------------------------------------------------------------------------- */
template<typename T>
inline void Math::matrixEig(__attribute__((unused)) UInt n,
__attribute__((unused)) T * A,
__attribute__((unused)) T * d,
__attribute__((unused)) T * V) {
AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function! Or implement it for the type " << debug::demangle(typeid(T).name()));
}
#ifdef AKANTU_USE_LAPACK
template<>
inline void Math::matrixEig<double>(UInt n, double * A, double * d, double * V) {
// Matrix A is row major, so the lapack function in fortran will process
// A^t. Asking for the left eigenvectors of A^t will give the transposed right
// eigenvectors of A so in the C++ code the right eigenvectors.
char jobvl;
if(V != NULL)
jobvl = 'V'; // compute left eigenvectors
else
jobvl = 'N'; // compute left eigenvectors
char jobvr('N'); // compute right eigenvectors
double * di = new double[n]; // imaginary part of the eigenvalues
int info;
int N = n;
double wkopt;
int lwork = -1;
// query and allocate the optimal workspace
dgeev_(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, NULL, &N, &wkopt, &lwork, &info);
lwork = int(wkopt);
double * work = new double[lwork];
// solve the eigenproblem
dgeev_(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, NULL, &N, work, &lwork, &info);
AKANTU_DEBUG_ASSERT(info == 0, "Problem computing eigenvalues/vectors. DGEEV exited with the value " << info);
delete [] work;
delete [] di; // I hope for you that there was no complex eigenvalues !!!
}
#endif
/* -------------------------------------------------------------------------- */
inline void Math::matrix22_eigenvalues(Real * A, Real *Adiag) {
///d = determinant of Matrix A
Real d = det2(A);
///b = trace of Matrix A
Real b = A[0]+A[3];
Real c = sqrt(b*b - 4 *d);
Adiag[0]= .5*(b + c);
Adiag[1]= .5*(b - c);
}
/* -------------------------------------------------------------------------- */
inline void Math::matrix33_eigenvalues(Real * A, Real *Adiag) {
/// a L^3 + b L^2 + c L + d = 0
matrixEig(3, A, Adiag);
// Real a = -1 ;
// ///b = trace of Matrix A
// Real b = A[0]+A[4]+A[8];
// /// c = 0.5*(trace(M^2)-trace(M)^2)
// Real c = A[1]*A[3] + A[2]*A[6] + A[5]*A[7] - A[0]*A[4] -
// A[0]*A[8] - A[4]*A[8];
// ///d = determinant of Matrix A
// Real d = det3(A);
//
// /// Define x, y, z
// Real x = c/a - b*b/(3.*a*a);
// Real y = 2.*b*b*b/(27.*a*a*a) - b*c/(3.*a*a) + d/a;
// Real z = y*y/4. + x*x*x/27.;
// /// Define I, j, k, m, n, p (so equations are not so cluttered)
// Real i = sqrt(y*y/4. - z);
// Real j = pow(i,1./3.);
// Real k = 0;
// if (std::abs(i) > 1e-12)
// k = acos(-(y/(2.*i)));
//
// Real m = cos(k/3);
// Real n = sqrt(3.)*sin(k/3);
// Real p = -b/(3.*a);
//
// Adiag[0] = 2*j*m + p;
// Adiag[1] = -j *(m + n) + p;
// Adiag[2] = -j * (m - n) + p;
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
inline void Math::eigenvalues(Real * A, Real * d) {
if(dim == 1) { d[0] = A[0]; }
else if(dim == 2) { matrix22_eigenvalues(A, d); }
// else if(dim == 3) { matrix33_eigenvalues(A, d); }
else matrixEig(dim, A, d);
}
/* -------------------------------------------------------------------------- */
inline Real Math::det2(const Real * mat) {
return mat[0]*mat[3] - mat[1]*mat[2];
}
/* -------------------------------------------------------------------------- */
inline Real Math::det3(const Real * mat) {
return
mat[0]*(mat[4]*mat[8]-mat[7]*mat[5])
- mat[3]*(mat[1]*mat[8]-mat[7]*mat[2])
+ mat[6]*(mat[1]*mat[5]-mat[4]*mat[2]);
}
/* -------------------------------------------------------------------------- */
template<UInt n>
inline Real Math::det(const Real * mat) {
if(n == 1) return *mat;
else if(n == 2) return det2(mat);
else if(n == 3) return det3(mat);
else return det(n, mat);
}
/* -------------------------------------------------------------------------- */
template<typename T>
inline T Math::det(__attribute__((unused)) UInt n,
__attribute__((unused)) const T * A) {
AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function!");
return T();
}
#ifdef AKANTU_USE_LAPACK
template<>
inline Real Math::det<double>(UInt n, const double * A) {
int N = n;
int info;
int * ipiv = new int[N+1];
double * LU = new double[N*N];
std::copy(A, A + N*N, LU);
// LU factorization of A
dgetrf_(&N, &N, LU, &N, ipiv, &info);
if(info > 0) {
AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: "
<< info <<" )");
}
// det(A) = det(L) * det(U) = 1 * det(U) = product_i U_{ii}
double det = 1.;
for (int i = 0; i < N; ++i) det *= (2*(ipiv[i] == i) - 1) * LU[i*n+i];
delete [] ipiv;
delete [] LU;
return det;
}
#endif
/* -------------------------------------------------------------------------- */
inline void Math::normal2(const Real * vec,Real * normal) {
normal[0] = vec[1];
normal[1] = -vec[0];
Math::normalize2(normal);
}
/* -------------------------------------------------------------------------- */
inline void Math::normal3(const Real * vec1,const Real * vec2,Real * normal) {
Math::vectorProduct3(vec1,vec2,normal);
Math::normalize3(normal);
}
/* -------------------------------------------------------------------------- */
inline void Math::normalize2(Real * vec) {
Real norm = Math::norm2(vec);
vec[0] /= norm;
vec[1] /= norm;
}
/* -------------------------------------------------------------------------- */
inline void Math::normalize3(Real * vec) {
Real norm = Math::norm3(vec);
vec[0] /= norm;
vec[1] /= norm;
vec[2] /= norm;
}
/* -------------------------------------------------------------------------- */
inline Real Math::norm2(const Real * vec) {
return sqrt(vec[0]*vec[0] + vec[1]*vec[1]);
}
/* -------------------------------------------------------------------------- */
inline Real Math::norm3(const Real * vec) {
return sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
}
/* -------------------------------------------------------------------------- */
inline Real Math::norm(UInt n, const Real * vec) {
Real norm = 0.;
for (UInt i = 0; i < n; ++i) {
norm += vec[i]*vec[i];
}
return sqrt(norm);
}
/* -------------------------------------------------------------------------- */
inline void Math::inv2(const Real * mat,Real * inv) {
Real det_mat = det2(mat);
inv[0] = mat[3] / det_mat;
inv[1] = -mat[1] / det_mat;
inv[2] = -mat[2] / det_mat;
inv[3] = mat[0] / det_mat;
}
/* -------------------------------------------------------------------------- */
inline void Math::inv3(const Real * mat,Real * inv) {
Real det_mat = det3(mat);
inv[0] = (mat[4]*mat[8] - mat[7]*mat[5])/det_mat;
inv[1] = (mat[2]*mat[7] - mat[8]*mat[1])/det_mat;
inv[2] = (mat[1]*mat[5] - mat[4]*mat[2])/det_mat;
inv[3] = (mat[5]*mat[6] - mat[8]*mat[3])/det_mat;
inv[4] = (mat[0]*mat[8] - mat[6]*mat[2])/det_mat;
inv[5] = (mat[2]*mat[3] - mat[5]*mat[0])/det_mat;
inv[6] = (mat[3]*mat[7] - mat[6]*mat[4])/det_mat;
inv[7] = (mat[1]*mat[6] - mat[7]*mat[0])/det_mat;
inv[8] = (mat[0]*mat[4] - mat[3]*mat[1])/det_mat;
}
/* -------------------------------------------------------------------------- */
template<UInt n>
inline void Math::inv(const Real * A, Real * Ainv) {
if(n == 1) *Ainv = 1./ *A;
else if(n == 2) inv2(A, Ainv);
else if(n == 3) inv3(A, Ainv);
else inv(n, A, Ainv);
}
/* -------------------------------------------------------------------------- */
template<typename T>
inline void Math::inv(__attribute__((unused)) UInt n,
__attribute__((unused)) const T * A,
__attribute__((unused)) T * Ainv) {
AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function! Or implement it for the type " << debug::demangle(typeid(T).name()));
}
#ifdef AKANTU_USE_LAPACK
template<>
inline void Math::inv<double>(UInt n, const double * A, double * invA) {
int N = n;
int info;
int * ipiv = new int[N+1];
int lwork = N*N;
double * work = new double[lwork];
std::copy(A, A + n*n, invA);
dgetrf_(&N, &N, invA, &N, ipiv, &info);
if(info > 0) {
AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: "
<< info <<" )");
}
dgetri_(&N, invA, &N, ipiv, work, &lwork, &info);
if(info != 0) {
AKANTU_DEBUG_ERROR("Cannot invert the matrix (info: "<< info <<" )");
}
delete [] ipiv;
delete [] work;
}
#endif
/* -------------------------------------------------------------------------- */
inline void Math::vectorProduct3(const Real * v1, const Real * v2, Real * res) {
res[0] = v1[1]*v2[2] - v1[2]*v2[1];
res[1] = v1[2]*v2[0] - v1[0]*v2[2];
res[2] = v1[0]*v2[1] - v1[1]*v2[0];
}
/* -------------------------------------------------------------------------- */
inline Real Math::vectorDot2(const Real * v1, const Real * v2) {
return (v1[0]*v2[0] + v1[1]*v2[1]);
}
/* -------------------------------------------------------------------------- */
inline Real Math::vectorDot3(const Real * v1, const Real * v2) {
return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]);
}
/* -------------------------------------------------------------------------- */
inline Real Math::distance_2d(const Real * x, const Real * y) {
return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1]));
}
/* -------------------------------------------------------------------------- */
inline Real Math::triangle_inradius(const Real * coord1,
const Real * coord2,
const Real * coord3) {
/**
* @f{eqnarray*}{
* r &=& A / s \\
* A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)} \\
* s &=& \frac{a + b + c}{2}
* @f}
*/
Real a, b, c;
a = distance_2d(coord1, coord2);
b = distance_2d(coord2, coord3);
c = distance_2d(coord1, coord3);
Real s;
s = (a + b + c) * 0.5;
return sqrt((s - a) * (s - b) * (s - c) / s);
}
/* -------------------------------------------------------------------------- */
inline Real Math::distance_3d(const Real * x, const Real * y) {
return sqrt((y[0] - x[0])*(y[0] - x[0])
+ (y[1] - x[1])*(y[1] - x[1])
+ (y[2] - x[2])*(y[2] - x[2])
);
}
/* -------------------------------------------------------------------------- */
inline Real Math::tetrahedron_volume(const Real * coord1,
const Real * coord2,
const Real * coord3,
const Real * coord4) {
Real xx[9], vol;
xx[0] = coord2[0]; xx[1] = coord2[1]; xx[2] = coord2[2];
xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2];
xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2];
vol = det3(xx);
xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2];
xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2];
xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2];
vol -= det3(xx);
xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2];
xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2];
xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2];
vol += det3(xx);
xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2];
xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2];
xx[6] = coord3[0]; xx[7] = coord3[1]; xx[8] = coord3[2];
vol -= det3(xx);
vol /= 6;
return vol;
}
/* -------------------------------------------------------------------------- */
inline Real Math::tetrahedron_inradius(const Real * coord1,
const Real * coord2,
const Real * coord3,
const Real * coord4) {
Real l12, l13, l14, l23, l24, l34;
l12 = distance_3d(coord1, coord2);
l13 = distance_3d(coord1, coord3);
l14 = distance_3d(coord1, coord4);
l23 = distance_3d(coord2, coord3);
l24 = distance_3d(coord2, coord4);
l34 = distance_3d(coord3, coord4);
Real s1, s2, s3, s4;
s1 = (l12 + l23 + l13) * 0.5;
s1 = sqrt(s1*(s1-l12)*(s1-l23)*(s1-l13));
s2 = (l12 + l24 + l14) * 0.5;
s2 = sqrt(s2*(s2-l12)*(s2-l24)*(s2-l14));
s3 = (l23 + l34 + l24) * 0.5;
s3 = sqrt(s3*(s3-l23)*(s3-l34)*(s3-l24));
s4 = (l13 + l34 + l14) * 0.5;
s4 = sqrt(s4*(s4-l13)*(s4-l34)*(s4-l14));
Real volume = Math::tetrahedron_volume(coord1,coord2,coord3,coord4);
return 3*volume/(s1+s2+s3+s4);
}
/* -------------------------------------------------------------------------- */
inline void Math::barycenter(const Real * coord,
UInt nb_points, UInt spatial_dimension,
Real * barycenter) {
memset(barycenter, 0, spatial_dimension * sizeof(Real));
for (UInt n = 0; n < nb_points; ++n) {
UInt offset = n * spatial_dimension;
for (UInt i = 0; i < spatial_dimension; ++i) {
barycenter[i] += coord[offset + i] / (Real) nb_points;
}
}
}
/* -------------------------------------------------------------------------- */
inline void Math::vector_2d(const Real * x, const Real * y, Real * res) {
res[0] = y[0]-x[0];
res[1] = y[1]-x[1];
}
/* -------------------------------------------------------------------------- */
inline void Math::vector_3d(const Real * x, const Real * y, Real * res) {
res[0] = y[0]-x[0];
res[1] = y[1]-x[1];
res[2] = y[2]-x[2];
}
/* -------------------------------------------------------------------------- */
inline bool Math::are_float_equal(const Real x, const Real y){
return (std::abs( x - y) < tolerance);
}
/* -------------------------------------------------------------------------- */
inline bool Math::isnan(Real x) {
#if defined(__INTEL_COMPILER)
#pragma warning ( push )
#pragma warning ( disable : 1572 )
#endif //defined(__INTEL_COMPILER)
// x = x return false means x = quiet_NaN
return !(x == x);
#if defined(__INTEL_COMPILER)
#pragma warning ( pop )
#endif //defined(__INTEL_COMPILER)
}
/* -------------------------------------------------------------------------- */
inline bool Math::are_vector_equal(UInt n, Real * x, Real * y){
bool test = true;
for (UInt i = 0; i < n; ++i) {
test &= are_float_equal(x[i],y[i]);
}
return test;
}
/* -------------------------------------------------------------------------- */
inline bool Math::intersects(Real x_min, Real x_max, Real y_min, Real y_max) {
- return ! ((x_max < y_min) || (x_min > y_max));
+ return ! ((x_max <= y_min) || (x_min >= y_max));
}
/* -------------------------------------------------------------------------- */
inline bool Math::is_in_range(Real a, Real x_min, Real x_max) {
return ((a >= x_min) && (a <= x_max));
}
diff --git a/src/fem/fem_template_cohesive.cc b/src/fem/fem_template_cohesive.cc
index 01ef5977c..9f1356416 100644
--- a/src/fem/fem_template_cohesive.cc
+++ b/src/fem/fem_template_cohesive.cc
@@ -1,169 +1,128 @@
/**
* @file fem_template_cohesive.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Wed Oct 31 16:24:42 2012
*
* @brief Specialization for cohesive element
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "fem_template.hh"
#include "shape_cohesive.hh"
#include "integrator_cohesive.hh"
__BEGIN_AKANTU__
-/* -------------------------------------------------------------------------- */
-// template <>
-// void FEMTemplate< IntegratorGauss, ShapeLagrange, _ek_cohesive >::
-// initShapeFunctions(const GhostType & ghost_type) {
-// AKANTU_DEBUG_IN();
-
-// UInt spatial_dimension = mesh->getSpatialDimension();
-// Mesh::type_iterator it = mesh->firstType(element_dimension, ghost_type, _ek_cohesive);
-// Mesh::type_iterator end = mesh->lastType(element_dimension, ghost_type, _ek_cohesive);
-// for(; it != end; ++it) {
-// ElementType type = *it;
-
-// #define INIT_SHAPE_FUNCTIONS(type) \
-// integrator.computeQuadraturePoints<type>(ghost_type); \
-// integrator. \
-// precomputeJacobiansOnQuadraturePoints<type>(ghost_type); \
-// integrator. \
-// checkJacobians<type>(ghost_type); \
-// const Vector<Real> & control_points = \
-// integrator.getQuadraturePoints<type>(ghost_type); \
-// shape_functions. \
-// setControlPointsByType<type>(control_points, ghost_type); \
-// shape_functions. \
-// precomputeShapesOnControlPoints<type>(ghost_type); \
-// if (element_dimension == spatial_dimension) \
-// shape_functions. \
-// precomputeShapeDerivativesOnControlPoints<type>(ghost_type);
-
-// AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS);
-// #undef INIT_SHAPE_FUNCTIONS
-// }
-// AKANTU_DEBUG_OUT();
-// }
-
-
/* -------------------------------------------------------------------------- */
/* compatibility functions */
/* -------------------------------------------------------------------------- */
template <>
Real FEMTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive>::integrate(const Vector<Real> & f,
const ElementType & type,
const GhostType & ghost_type,
const Vector<UInt> * filter_elements) const{
AKANTU_DEBUG_IN();
#ifndef AKANTU_NDEBUG
-// std::stringstream sstr; sstr << ghost_type;
-// AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
-// "The vector " << nablauq.getID() << " is not taged " << ghost_type);
UInt nb_element = mesh->getNbElement(type, ghost_type);
if(filter_elements != NULL) nb_element = filter_elements->getSize();
UInt nb_quadrature_points = getNbQuadraturePoints(type);
AKANTU_DEBUG_ASSERT(f.getSize() == nb_element * nb_quadrature_points,
"The vector f(" << f.getID()
<< ") has not the good size.");
AKANTU_DEBUG_ASSERT(f.getNbComponent() == 1,
"The vector f(" << f.getID()
<< ") has not the good number of component.");
#endif
Real integral = 0.;
#define INTEGRATE(type) \
integral = integrator. integrate<type>(f, \
ghost_type, \
filter_elements);
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INTEGRATE);
#undef INTEGRATE
AKANTU_DEBUG_OUT();
return integral;
}
/* -------------------------------------------------------------------------- */
template <>
void FEMTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive>
::integrate(const Vector<Real> & f,
Vector<Real> &intf,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type,
const Vector<UInt> * filter_elements) const{
#ifndef AKANTU_NDEBUG
-// std::stringstream sstr; sstr << ghost_type;
-// AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
-// "The vector " << nablauq.getID() << " is not taged " << ghost_type);
UInt nb_element = mesh->getNbElement(type, ghost_type);
if(filter_elements != NULL) nb_element = filter_elements->getSize();
UInt nb_quadrature_points = getNbQuadraturePoints(type);
AKANTU_DEBUG_ASSERT(f.getSize() == nb_element * nb_quadrature_points,
"The vector f(" << f.getID() << " size " << f.getSize()
<< ") has not the good size (" << nb_element << ").");
AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom ,
"The vector f(" << f.getID()
<< ") has not the good number of component.");
AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
"The vector intf(" << intf.getID()
<< ") has not the good number of component.");
AKANTU_DEBUG_ASSERT(intf.getSize() == nb_element,
"The vector intf(" << intf.getID()
<< ") has not the good size.");
#endif
#define INTEGRATE(type) \
integrator. integrate<type>(f, \
intf, \
nb_degree_of_freedom, \
ghost_type, \
filter_elements);
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INTEGRATE);
#undef INTEGRATE
}
/* -------------------------------------------------------------------------- */
template <>
void FEMTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive>::
gradientOnQuadraturePoints(__attribute__((unused)) const Vector<Real> &u,
__attribute__((unused)) Vector<Real> &nablauq,
__attribute__((unused)) const UInt nb_degree_of_freedom,
__attribute__((unused)) const ElementType & type,
__attribute__((unused)) const GhostType & ghost_type,
__attribute__((unused)) const Vector<UInt> * filter_elements) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
__END_AKANTU__
diff --git a/src/fem/geometrical_data_tmpl.hh b/src/fem/geometrical_data_tmpl.hh
index 43ee844f3..eca6c586a 100644
--- a/src/fem/geometrical_data_tmpl.hh
+++ b/src/fem/geometrical_data_tmpl.hh
@@ -1,194 +1,194 @@
/**
* @file geometrical_data_tmpl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Wed Nov 14 14:57:27 2012
*
* @brief Specialization of the geometrical types
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_point> {
static const UInt spatial_dimension = 0;
static const UInt nb_nodes_per_element = 1;
static const GeometricalType p1_element_type = _gt_point;
static const GeometricalType facet_type = _gt_not_defined;
static const UInt nb_facets = 0;
static const UInt * facet_connectivity[] = {};
};
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_segment_2> {
static const UInt spatial_dimension = 1;
static const UInt nb_nodes_per_element = 2;
static const GeometricalType p1_element_type = _gt_segment_2;
static const UInt nb_facets = 2;
static const GeometricalType facet_type = _gt_point;
static const UInt vec_facet_connectivity[] = {0,
- 1};
+ 1};
static const UInt * facet_connectivity[] = {&vec_facet_connectivity[0],
- &vec_facet_connectivity[1]};
+ &vec_facet_connectivity[1]};
};
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_segment_3> {
static const UInt spatial_dimension = 1;
static const UInt nb_nodes_per_element = 3;
static const GeometricalType p1_element_type = _gt_segment_2;
static const UInt nb_facets = 2;
static const GeometricalType facet_type = _gt_point;
static const UInt vec_facet_connectivity[] = {0,
- 1};
+ 1};
static const UInt * facet_connectivity[] = {&vec_facet_connectivity[0],
- &vec_facet_connectivity[1]};
+ &vec_facet_connectivity[1]};
};
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_triangle_3> {
static const UInt spatial_dimension = 2;
static const UInt nb_nodes_per_element = 3;
static const GeometricalType p1_element_type = _gt_triangle_3;
static const UInt nb_facets = 3;
static const GeometricalType facet_type = _gt_segment_2;
static const UInt vec_facet_connectivity[] = {0, 1,
- 1, 2,
- 2, 0};
+ 1, 2,
+ 2, 0};
static const UInt * facet_connectivity[] = {&vec_facet_connectivity[0],
- &vec_facet_connectivity[2],
- &vec_facet_connectivity[4]};
+ &vec_facet_connectivity[2],
+ &vec_facet_connectivity[4]};
};
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_triangle_6> {
static const UInt spatial_dimension = 2;
static const UInt nb_nodes_per_element = 6;
static const GeometricalType p1_element_type = _gt_triangle_3;
static const UInt nb_facets = 3;
static const GeometricalType facet_type = _gt_segment_3;
static const UInt vec_facet_connectivity[] = {0, 1, 3,
- 1, 2, 4,
- 2, 0, 5};
+ 1, 2, 4,
+ 2, 0, 5};
static const UInt * facet_connectivity[] = {&vec_facet_connectivity[0],
- &vec_facet_connectivity[3],
- &vec_facet_connectivity[6]};
+ &vec_facet_connectivity[3],
+ &vec_facet_connectivity[6]};
};
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_tetrahedron_4> {
static const UInt spatial_dimension = 3;
static const UInt nb_nodes_per_element = 4;
static const GeometricalType p1_element_type = _gt_tetrahedron_4;
static const UInt nb_facets = 4;
static const GeometricalType facet_type = _gt_triangle_3;
static const UInt vec_facet_connectivity[] = {0, 2, 1,
- 1, 2, 3,
- 2, 0, 3,
- 0, 1, 3};
+ 1, 2, 3,
+ 2, 0, 3,
+ 0, 1, 3};
static const UInt * facet_connectivity[] = {&vec_facet_connectivity[0],
- &vec_facet_connectivity[3],
- &vec_facet_connectivity[6],
- &vec_facet_connectivity[9]};
+ &vec_facet_connectivity[3],
+ &vec_facet_connectivity[6],
+ &vec_facet_connectivity[9]};
};
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_tetrahedron_10> {
static const UInt spatial_dimension = 3;
static const UInt nb_nodes_per_element = 10;
static const GeometricalType p1_element_type = _gt_tetrahedron_4;
static const UInt nb_facets = 4;
static const GeometricalType facet_type = _gt_triangle_6;
static const UInt vec_facet_connectivity[] = {0, 2, 1, 6, 5, 4,
- 1, 2, 3, 5, 9, 8,
- 2, 0, 3, 6, 7, 9,
- 0, 1, 3, 4, 8, 7};
+ 1, 2, 3, 5, 9, 8,
+ 2, 0, 3, 6, 7, 9,
+ 0, 1, 3, 4, 8, 7};
static const UInt * facet_connectivity[] = {&vec_facet_connectivity[0],
- &vec_facet_connectivity[6],
- &vec_facet_connectivity[12],
- &vec_facet_connectivity[18]};
+ &vec_facet_connectivity[6],
+ &vec_facet_connectivity[12],
+ &vec_facet_connectivity[18]};
};
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_quadrangle_4> {
static const UInt spatial_dimension = 2;
static const UInt nb_nodes_per_element = 4;
static const GeometricalType p1_element_type = _gt_quadrangle_4;
static const UInt nb_facets = 4;
static const GeometricalType facet_type = _gt_segment_2;
static const UInt vec_facet_connectivity[] = {0, 1,
- 1, 2,
- 2, 3,
- 3, 0};
+ 1, 2,
+ 2, 3,
+ 3, 0};
static const UInt * facet_connectivity[] = {&vec_facet_connectivity[0],
- &vec_facet_connectivity[2],
- &vec_facet_connectivity[4],
- &vec_facet_connectivity[6]};
+ &vec_facet_connectivity[2],
+ &vec_facet_connectivity[4],
+ &vec_facet_connectivity[6]};
};
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_quadrangle_8> {
static const UInt spatial_dimension = 2;
static const UInt nb_nodes_per_element = 8;
static const GeometricalType p1_element_type = _gt_quadrangle_4;
static const UInt nb_facets = 4;
static const GeometricalType facet_type = _gt_segment_3;
static const UInt vec_facet_connectivity[] = {0, 1, 4,
- 1, 2, 5,
- 2, 3, 6,
- 3, 0, 7};
+ 1, 2, 5,
+ 2, 3, 6,
+ 3, 0, 7};
static const UInt * facet_connectivity[] = {vec_facet_connectivity + 0,
- vec_facet_connectivity + 3,
- vec_facet_connectivity + 6,
- vec_facet_connectivity + 9};
+ vec_facet_connectivity + 3,
+ vec_facet_connectivity + 6,
+ vec_facet_connectivity + 9};
};
/* -------------------------------------------------------------------------- */
template<>
struct GeometricalElementData<_gt_hexahedron_8> {
static const UInt spatial_dimension = 3;
static const UInt nb_nodes_per_element = 8;
static const GeometricalType p1_element_type = _gt_hexahedron_8;
static const UInt nb_facets = 6;
static const GeometricalType facet_type = _gt_quadrangle_4;
static const UInt vec_facet_connectivity[] = {0, 1, 2, 3,
- 0, 1, 5, 4,
- 1, 2, 6, 5,
- 2, 3, 7, 6,
- 3, 0, 4, 7,
- 4, 5, 6, 7};
+ 0, 1, 5, 4,
+ 1, 2, 6, 5,
+ 2, 3, 7, 6,
+ 3, 0, 4, 7,
+ 4, 5, 6, 7};
static const UInt * facet_connectivity[] = {&vec_facet_connectivity[0],
- &vec_facet_connectivity[4],
- &vec_facet_connectivity[8],
- &vec_facet_connectivity[12],
- &vec_facet_connectivity[16],
- &vec_facet_connectivity[20]};
+ &vec_facet_connectivity[4],
+ &vec_facet_connectivity[8],
+ &vec_facet_connectivity[12],
+ &vec_facet_connectivity[16],
+ &vec_facet_connectivity[20]};
};
diff --git a/src/io/dumper/dumper_iohelper.hh b/src/io/dumper/dumper_iohelper.hh
index d7a3ec310..387613686 100644
--- a/src/io/dumper/dumper_iohelper.hh
+++ b/src/io/dumper/dumper_iohelper.hh
@@ -1,207 +1,216 @@
/**
* @file dumper_iohelper.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Fri Oct 26 21:52:40 2012
*
* @brief Define the akantu dumper interface for IOhelper dumpers
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_DUMPER_IOHELPER_HH__
#define __AKANTU_DUMPER_IOHELPER_HH__
#include "aka_common.hh"
#include "aka_types.hh"
#include "mesh.hh"
#include <io_helper.hh>
__BEGIN_AKANTU__
class DumperIOHelper {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
DumperIOHelper();
virtual ~DumperIOHelper();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
class Field;
void registerMesh(const Mesh & mesh, UInt spatial_dimension = 0,
const GhostType & ghost_type = _not_ghost,
const ElementKind & element_kind = _ek_not_defined);
void registerField(const std::string & field_id, Field * field);
void unRegisterField(const std::string & field_id);
void dump();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(Dumper, *dumper, iohelper::Dumper &)
static iohelper::ElemType getIOHelperType(ElementType type);
protected:
template <ElementType type> static iohelper::ElemType getIOHelperType();
public:
/* ------------------------------------------------------------------------ */
/* Field descriptors */
/* ------------------------------------------------------------------------ */
/// Field interface
class Field {
public:
Field() : padding_n(0), padding_m(0) {}
virtual ~Field() {};
virtual void registerToDumper(const std::string & id, iohelper::Dumper & dupmer) = 0;
virtual void setPadding(UInt n, UInt m = 1) {
padding_n = n;
padding_m = m;
}
protected:
UInt padding_n, padding_m;
};
/* ------------------------------------------------------------------------ */
template<class T, template<class> class R>
class iterator_helper;
template<class T, template<class> class R>
class PaddingHelper;
/* ------------------------------------------------------------------------ */
/* Nodal field wrapper */
template<typename T>
class NodalField;
/* ------------------------------------------------------------------------ */
/* Generic class used as interface for the others */
template<typename i_type, typename d_type,
template<typename> class ret_type, class daughter>
class element_iterator;
template<typename i_type, typename d_type,
template<typename> class ret_type, class daughter>
+ class generic_quadrature_point_iterator;
+
+ template<typename T, template<typename> class ret_type>
class quadrature_point_iterator;
template<typename T, class iterator_type,
template<typename> class ret_type>
class GenericElementalField;
- template<typename T, class iterator_type,
+ template<typename T,
+ class iterator_type,
template<typename> class ret_type>
class GenericQuadraturePointsField;
+ template<typename T,
+ template<typename> class ret_type,
+ template<typename, template<class> class> class iterator_type = quadrature_point_iterator>
+ class QuadraturePointsField;
+
/* ------------------------------------------------------------------------ */
/* Elemental Fields wrapper */
class element_type_field_iterator;
class element_partition_field_iterator;
template<typename T, template<class> class ret_type>
class elemental_field_iterator;
class ElementTypeField;
class ElementPartitionField;
template<typename T, template<typename> class ret_type = types::Vector>
class ElementalField;
/* ------------------------------------------------------------------------ */
/* Material Field wrapper */
template<typename T,
template<class> class ret_type,
template<typename, template<class> class> class padding_helper_type,
template<typename, template<class> class> class int_iterator>
class generic_internal_material_field_iterator;
template<typename T, template<class> class ret_type>
class internal_material_field_iterator;
template<typename T, template<class> class ret_type>
class material_stress_field_iterator;
template<typename T, template<class> class ret_type>
class material_strain_field_iterator;
template<typename T, template<class> class ret_type = types::Vector,
template<typename, template<class> class> class iterator_type = internal_material_field_iterator>
class InternalMaterialField;
template<class T, template<class> class R>
class MaterialPaddingHelper;
template<class T, template<class> class R>
class StressPaddingHelper;
template<class T, template<class> class R>
class StrainPaddingHelper;
/* ------------------------------------------------------------------------ */
/* Field homogenizing wrapper */
template<typename T, class Container, template<class> class sub_type>
class PaddingHomogenizingFunctor;
template<typename T, class Container, template<class> class sub_type>
class AvgHomogenizingFunctor;
template<typename T, template< typename,
template<class> class,
template<typename, template<class> class> class > class Container,
+ template<typename, template<class> class> class sub_iterator = internal_material_field_iterator,
template<typename, class, template<class> class> class Funct = AvgHomogenizingFunctor,
- template<typename> class ret_type = types::Vector,
- template<typename, template<class> class> class sub_iterator = internal_material_field_iterator>
+ template<typename> class ret_type = types::Vector>
class HomogenizedField;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// internal iohelper::Dumper
iohelper::Dumper * dumper;
typedef std::map<std::string, Field *> Fields;
/// list of registered fields to dump
Fields fields;
/// dump counter
UInt count;
/// filename
std::string filename;
};
#include "dumper_iohelper_tmpl.hh"
__END_AKANTU__
#endif /* __AKANTU_DUMPER_IOHELPER_HH__ */
diff --git a/src/io/dumper/dumper_iohelper_tmpl.hh b/src/io/dumper/dumper_iohelper_tmpl.hh
index 8c22819db..a10acf47e 100644
--- a/src/io/dumper/dumper_iohelper_tmpl.hh
+++ b/src/io/dumper/dumper_iohelper_tmpl.hh
@@ -1,96 +1,94 @@
/**
* @file dumper_iohelper_tmpl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Fri Oct 26 21:52:40 2012
*
* @brief implementation of the DumperIOHelper class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template<class T, template<class> class R>
class DumperIOHelper::iterator_helper {
typedef typename Vector<T>::template const_iterator< R<T> > internal_iterator;
public:
static internal_iterator begin(const Vector<T> & vect, UInt m, UInt n, UInt size) {
return vect.begin_reinterpret(n*m, size);
}
static internal_iterator end(const Vector<T> & vect, UInt m, UInt n, UInt size) {
return vect.end_reinterpret(n*m, size);
}
};
template<class T>
class DumperIOHelper::iterator_helper<T, types::Matrix> {
typedef typename Vector<T>::template const_iterator< types::Matrix<T> > internal_iterator;
public:
static internal_iterator begin(const Vector<T> & vect, UInt m, UInt n, UInt size) {
return vect.begin_reinterpret(m, n, size);
}
static internal_iterator end(const Vector<T> & vect, UInt m, UInt n, UInt size) {
return vect.end_reinterpret(m, n, size);
}
};
/* -------------------------------------------------------------------------- */
template<class T, template<class> class R>
class DumperIOHelper::PaddingHelper {
public:
- inline R<T> pad(const R<T> & in,
+ static inline R<T> pad(const R<T> & in,
__attribute__((unused)) UInt padding_m,
__attribute__((unused)) UInt padding_n,
__attribute__((unused)) UInt nb_data) {
return in; // trick due to the fact that IOHelper padds the vectors (avoid a copy of data)
}
};
template<class T>
class DumperIOHelper::PaddingHelper<T, types::Matrix> {
public:
- inline types::Matrix<T> pad(const types::Matrix<T> & in, UInt padding_m, UInt padding_n, UInt nb_data) {
+ static inline types::Matrix<T> pad(const types::Matrix<T> & in, UInt padding_m, UInt padding_n, UInt nb_data) {
if(padding_m <= in.rows() && padding_n*nb_data <= in.cols())
return in;
else {
types::Matrix<T> ret(padding_m, padding_n * nb_data);
for (UInt d = 0; d < nb_data; ++d)
for (UInt i = 0; i < in.rows(); ++i)
for (UInt j = 0; j < in.cols() / nb_data; ++j)
ret(i, j + d * padding_n) = in(i, j + d * in.cols());
return ret;
}
}
};
#include "dumper_iohelper_tmpl_nodal_field.hh"
-
#include "dumper_iohelper_tmpl_elemental_field.hh"
#include "dumper_iohelper_tmpl_quadrature_points_field.hh"
-
diff --git a/src/io/dumper/dumper_iohelper_tmpl_homogenizing_field.hh b/src/io/dumper/dumper_iohelper_tmpl_homogenizing_field.hh
index 5da1553b5..28235fd80 100644
--- a/src/io/dumper/dumper_iohelper_tmpl_homogenizing_field.hh
+++ b/src/io/dumper/dumper_iohelper_tmpl_homogenizing_field.hh
@@ -1,224 +1,249 @@
/**
* @file dumper_iohelper_tmpl_homogenizing_field.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Fri Oct 26 21:52:40 2012
*
* @brief description of field homogenizing field
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
__END_AKANTU__
#include "solid_mechanics_model.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
/* Homogenizing containers */
/* -------------------------------------------------------------------------- */
template<typename T, class Container, template<class> class sub_type>
class DumperIOHelper::PaddingHomogenizingFunctor {
public:
PaddingHomogenizingFunctor(Container & cont) : cont(cont){
}
UInt getNbComponent() {
typename Container::iterator it = cont.begin();
typename Container::iterator end = cont.end();
UInt nb_comp = 0;
for (; it != end; ++it) nb_comp = std::max(nb_comp, (*it).size());
return nb_comp;
}
sub_type<T> operator()(const sub_type<T> & vect, __attribute__((unused)) const ElementType & type) {
return vect;
}
protected:
Container & cont;
};
/* -------------------------------------------------------------------------- */
template<typename T, class Container, template<class> class sub_type>
class DumperIOHelper::AvgHomogenizingFunctor {
public:
AvgHomogenizingFunctor(Container & cont) : cont(cont){
}
UInt getNbComponent() {
typename Container::iterator it = cont.begin();
typename Container::iterator end = cont.end();
UInt nb_comp = 0;
UInt i = 0;
for (; it != end; ++it) {
ElementType type = it.getType();
UInt nb_quad = cont.getFEM().getNbQuadraturePoints(type);
nb_comp = std::max(nb_comp, (*it).size() / nb_quad);
++i;
}
return nb_comp;
}
sub_type<T> operator()(const sub_type<T> & vect, const ElementType & type) {
UInt nb_quad = cont.getFEM().getNbQuadraturePoints(type);
if(nb_quad == 1) {
return vect;
} else {
UInt s = vect.size() / nb_quad;
sub_type<T> v(s);
for (UInt i = 0; i < s; ++i) {
for (UInt q = 0; q < nb_quad; ++q) {
v[i] += vect[i + q*s] / Real(nb_quad);
}
}
return v;
}
}
protected:
Container & cont;
};
/* specialization */
template<typename T, class Container>
class DumperIOHelper::AvgHomogenizingFunctor<T, Container, types::Matrix> {
public:
AvgHomogenizingFunctor(Container & cont) : cont(cont){
}
UInt getNbComponent() {
typename Container::iterator it = cont.begin();
typename Container::iterator end = cont.end();
UInt nb_comp = 0;
UInt i = 0;
for (; it != end; ++it) {
ElementType type = it.getType();
UInt nb_quad = cont.getFEM().getNbQuadraturePoints(type);
nb_comp = std::max(nb_comp, (*it).size() / nb_quad);
++i;
}
return nb_comp;
}
types::Matrix<T> operator()(const types::Matrix<T> & vect,
const ElementType & type) {
UInt nb_quad = cont.getFEM().getNbQuadraturePoints(type);
if(nb_quad == 1) {
return vect;
} else {
UInt m = vect.rows();
UInt n = vect.cols() / nb_quad;
types::Matrix<T> v(m, n);
for (UInt q = 0; q < nb_quad; ++q)
for (UInt i = 0; i < m; ++i)
for (UInt j = 0; j < n; ++j)
v(i, j) += vect(i, j + q*n) / Real(nb_quad);
return v;
}
}
protected:
Container & cont;
};
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
-template<typename T, template< typename, template<class> class, template<typename, template<class> class> class> class Container,
+template<typename T,
+ template< typename,
+ template<class> class,
+ template<typename, template<class> class> class> class Container,
+ template<typename, template<class> class> class int_sub_iterator,
template<typename, class, template<class> class> class Funct,
- template<typename> class ret_type,
- template<typename, template<class> class> class int_sub_iterator>
+ template<typename> class ret_type>
class DumperIOHelper::HomogenizedField : public Field {
protected:
typedef typename Container<T, ret_type, int_sub_iterator>::iterator sub_iterator;
typedef Funct<T, Container<T, ret_type, int_sub_iterator>, ret_type> Functor;
public:
class iterator {
public:
iterator(const sub_iterator & it, Functor & func) : it(it), func(func) {}
bool operator!=(const iterator & it) { return it.it != this->it; }
iterator operator++() { ++this->it; return *this; }
ret_type<T> operator*() {
return func(*it, it.getType());
}
protected:
sub_iterator it;
Functor & func;
};
HomogenizedField(const SolidMechanicsModel & model,
const std::string & field_id,
UInt spatial_dimension = 0,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) :
cont(model, field_id, spatial_dimension, ghost_type, element_kind),
funct(cont) {
nb_component = funct.getNbComponent();
}
HomogenizedField(const SolidMechanicsModel & model,
const std::string & field_id,
UInt n,
UInt spatial_dimension = 0,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) :
cont(model, field_id, n, spatial_dimension, ghost_type, element_kind),
funct(cont) {
nb_component = funct.getNbComponent();
}
+ HomogenizedField(const FEM & fem,
+ const ByElementTypeVector<T> & field,
+ UInt spatial_dimension = 0,
+ GhostType ghost_type = _not_ghost,
+ ElementKind element_kind = _ek_not_defined) :
+ cont(fem, field, spatial_dimension, ghost_type, element_kind),
+ funct(cont) {
+ nb_component = funct.getNbComponent();
+ }
+
+ HomogenizedField(const FEM & fem,
+ const ByElementTypeVector<T> & field,
+ UInt n,
+ UInt spatial_dimension = 0,
+ GhostType ghost_type = _not_ghost,
+ ElementKind element_kind = _ek_not_defined) :
+ cont(fem, field, n, spatial_dimension, ghost_type, element_kind),
+ funct(cont) {
+ nb_component = funct.getNbComponent();
+ }
+
+
iterator begin() { return iterator(cont.begin(), funct); }
iterator end () { return iterator(cont.end(), funct); }
virtual void registerToDumper(const std::string & id,
iohelper::Dumper & dupmer) {
dupmer.addElemDataField(id, *this);
}
UInt getDim() {
if(padding_n && padding_m)
return padding_m*padding_n;
else return nb_component;
}
UInt size() { return cont.size(); }
iohelper::DataType getDataType() { return iohelper::getDataType<T>(); }
UInt isHomogeneous() { return true; }
void setPadding(UInt n, UInt m = 1) {
Field::setPadding(n, m);
cont.setPadding(n, m);
}
protected:
Container<T, ret_type, int_sub_iterator> cont;
Functor funct;
UInt nb_component;
};
diff --git a/src/io/dumper/dumper_iohelper_tmpl_material_internal_field.hh b/src/io/dumper/dumper_iohelper_tmpl_material_internal_field.hh
index 56d009657..9b125c07f 100644
--- a/src/io/dumper/dumper_iohelper_tmpl_material_internal_field.hh
+++ b/src/io/dumper/dumper_iohelper_tmpl_material_internal_field.hh
@@ -1,375 +1,375 @@
/**
* @file dumper_iohelper_tmpl_material_internal_field.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Fri Oct 26 21:52:40 2012
*
* @brief description of material internal field
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template<class T, template<class> class R>
class DumperIOHelper::MaterialPaddingHelper : public PaddingHelper<T, R> {
public:
MaterialPaddingHelper(const SolidMechanicsModel & model) : model(model) { }
inline R<T> pad(const R<T> & in,
UInt padding_m, UInt padding_n,
UInt nb_data,
__attribute__((unused)) UInt material_id) {
return PaddingHelper<T, R>::pad(in, padding_m, padding_n, nb_data);
}
protected:
const SolidMechanicsModel & model;
};
/* -------------------------------------------------------------------------- */
template<class T, template<class> class R>
class DumperIOHelper::StressPaddingHelper : public MaterialPaddingHelper<T, R> {
public:
StressPaddingHelper(const SolidMechanicsModel & model) : MaterialPaddingHelper<T, R>(model) {
AKANTU_DEBUG_ERROR("This class exists only to pad stress (types::Matrix<Real>) to 3D");
}
inline R<T> pad(const R<T> & in, UInt padding_m, UInt padding_n, UInt nb_data, UInt material_id) {}
protected:
const Material & material;
};
template<>
class DumperIOHelper::StressPaddingHelper<Real,
types::Matrix> : public MaterialPaddingHelper<Real,
types::Matrix> {
public:
StressPaddingHelper(const SolidMechanicsModel & model) :
MaterialPaddingHelper<Real, types::Matrix>(model) { }
inline types::Matrix<Real> pad(const types::Matrix<Real> & in,
UInt padding_m, UInt padding_n, UInt nb_data,
UInt material_id) {
if(padding_m <= in.rows() && padding_n * nb_data <= in.cols())
return in;
else {
AKANTU_DEBUG_ASSERT(padding_m == 3 && padding_n == 3, "This function can only pad to 3D");
types::Matrix<Real> stress =
MaterialPaddingHelper<Real, types::Matrix>::pad(in, 3, 3, nb_data, material_id);
if(in.rows() == 2 && in.cols() == 2 * nb_data) {
const Material & material = model.getMaterial(material_id);
bool plane_strain = !material.getParam<bool>("Plane_Stress");
if(plane_strain) {
Real nu = material.getParam<Real>("nu");
for (UInt d = 0; d < nb_data; ++d) {
stress(2, 2 + 3*d) = nu * (stress(0, 0 + 3*d) + stress(1, 1 + 3*d));
}
}
}
return stress;
}
}
};
/* -------------------------------------------------------------------------- */
template<class T, template<class> class R>
class DumperIOHelper::StrainPaddingHelper : public MaterialPaddingHelper<T, R> {
public:
StrainPaddingHelper(const SolidMechanicsModel & model) : MaterialPaddingHelper<T, R>(model) {
AKANTU_DEBUG_ERROR("This class exists only to pad strain (types::Matrix<Real>) to 3D");
}
inline R<T> pad(const R<T> & in, UInt padding_m, UInt padding_n, UInt nb_data, UInt material_id) {}
};
template<>
class DumperIOHelper::StrainPaddingHelper<Real, types::Matrix> : public MaterialPaddingHelper<Real, types::Matrix> {
public:
StrainPaddingHelper(const SolidMechanicsModel & model) :
MaterialPaddingHelper<Real, types::Matrix>(model) { }
inline types::Matrix<Real> pad(const types::Matrix<Real> & in,
UInt padding_m, UInt padding_n, UInt nb_data,
UInt material_id) {
if(padding_m <= in.rows() && padding_n * nb_data <= in.cols())
return in;
else {
AKANTU_DEBUG_ASSERT(padding_m == 3 && padding_n == 3, "This function can only pad to 3D");
types::Matrix<Real> strain =
MaterialPaddingHelper<Real, types::Matrix>::pad(in, 3, 3, nb_data, material_id);
if(in.rows() == 2 && in.cols() == 2 * nb_data) {
const Material & material = model.getMaterial(material_id);
bool plane_stress = material.getParam<bool>("Plane_Stress");
if(plane_stress) {
Real nu = material.getParam<Real>("nu");
for (UInt d = 0; d < nb_data; ++d) {
strain(2, 2 + 3*d) = nu / (nu - 1) * (strain(0, 0 + 3*d) + strain(1, 1 + 3*d));
}
}
}
return strain;
}
}
};
/* -------------------------------------------------------------------------- */
/* Element material field iterator/container */
/* -------------------------------------------------------------------------- */
template<typename T,
template<class> class ret_type,
template<typename, template<class> class> class padding_helper_type,
template<typename, template<class> class> class int_iterator>
-class DumperIOHelper::generic_internal_material_field_iterator : public quadrature_point_iterator< UInt, T, ret_type,
- int_iterator<T, ret_type> > {
+class DumperIOHelper::generic_internal_material_field_iterator : public generic_quadrature_point_iterator< UInt, T, ret_type,
+ int_iterator<T, ret_type> > {
public:
- typedef quadrature_point_iterator<UInt, T, ret_type, int_iterator<T, ret_type> > parent;
+ typedef generic_quadrature_point_iterator<UInt, T, ret_type, int_iterator<T, ret_type> > parent;
typedef typename parent::it_type it_type;
typedef typename parent::data_type data_type;
typedef typename parent::return_type return_type;
typedef typename parent::field_type field_type;
typedef typename parent::internal_iterator internal_iterator;
typedef typename Vector<T>::template const_iterator< ret_type<T> > internal_material_iterator;
public:
generic_internal_material_field_iterator(const field_type & element_material,
UInt n,
const typename field_type::type_iterator & t_it,
const typename field_type::type_iterator & t_it_end,
const internal_iterator & it,
ElementType element_type,
const GhostType ghost_type = _not_ghost) :
parent(element_material, 2, t_it, t_it_end,
it, element_type, ghost_type),
out_n(n), model(NULL), padding_helper(NULL) { }
~generic_internal_material_field_iterator() { delete padding_helper; }
generic_internal_material_field_iterator(const generic_internal_material_field_iterator & it) : parent(it) {
out_n = it.out_n;
field_id = it.field_id;
if(it.model) {
model = it.model;
padding_helper = new padding_helper_type<T, ret_type>(*model);
}
}
return_type operator*() {
const ret_type<UInt> & material_id = *this->vit;
UInt nb_data = this->fem->getNbQuadraturePoints(*this->tit);
try {
const Vector<T> & vect =
model->getMaterial(material_id[1]).getVector(field_id,
*this->tit,
this->ghost_type);
if(vect.getSize() == 0 || vect.getSize() < material_id[0]) // vector exists but has a wrong size
return return_type();
UInt ln = out_n;
if(out_n == 0) ln = vect.getNbComponent();
internal_material_iterator it
= iterator_helper<T, ret_type>::begin(vect, ln,
vect.getNbComponent() / ln * nb_data,
vect.getSize() / nb_data);
it += material_id[0];
return padding_helper->pad(*it, this->padding_m, this->padding_n, nb_data, material_id[1]);
} catch (...) {
return return_type();
}
}
void setModel(const SolidMechanicsModel & model) {
this->model = &model;
this->padding_helper = new padding_helper_type<T, ret_type>(model);
}
void setFieldID(const std::string & field_id) { this->field_id = field_id; }
protected:
virtual UInt getNbDataPerElem(__attribute__((unused)) const ElementType & type) { return 1; }
protected:
UInt out_n;
const SolidMechanicsModel * model;
std::string field_id;
padding_helper_type<T, ret_type> * padding_helper;
};
/* -------------------------------------------------------------------------- */
template<typename T, template<class> class ret_type>
class DumperIOHelper::internal_material_field_iterator :
public generic_internal_material_field_iterator<T, ret_type,
MaterialPaddingHelper,
internal_material_field_iterator> {
public:
typedef generic_internal_material_field_iterator<T, ret_type,
MaterialPaddingHelper,
DumperIOHelper::internal_material_field_iterator> parent;
typedef typename parent::it_type it_type;
typedef typename parent::data_type data_type;
typedef typename parent::return_type return_type;
typedef typename parent::field_type field_type;
typedef typename parent::internal_iterator internal_iterator;
typedef typename Vector<T>::template const_iterator< ret_type<T> > internal_material_iterator;
public:
internal_material_field_iterator(const field_type & element_material,
UInt n,
const typename field_type::type_iterator & t_it,
const typename field_type::type_iterator & t_it_end,
const internal_iterator & it,
ElementType element_type,
const GhostType ghost_type = _not_ghost) :
parent(element_material, n, t_it, t_it_end,
it, element_type, ghost_type) { }
};
/* -------------------------------------------------------------------------- */
template<typename T, template<class> class ret_type>
class DumperIOHelper::material_stress_field_iterator :
public generic_internal_material_field_iterator<T, ret_type, StressPaddingHelper,
material_stress_field_iterator> {
public:
typedef generic_internal_material_field_iterator<T, ret_type, StressPaddingHelper,
DumperIOHelper::material_stress_field_iterator> parent;
typedef typename parent::it_type it_type;
typedef typename parent::data_type data_type;
typedef typename parent::return_type return_type;
typedef typename parent::field_type field_type;
typedef typename parent::internal_iterator internal_iterator;
typedef typename Vector<T>::template const_iterator< ret_type<T> > internal_material_iterator;
public:
material_stress_field_iterator(const field_type & element_material,
UInt n,
const typename field_type::type_iterator & t_it,
const typename field_type::type_iterator & t_it_end,
const internal_iterator & it,
ElementType element_type,
const GhostType ghost_type = _not_ghost) :
parent(element_material, n, t_it, t_it_end,
it, element_type, ghost_type) { }
};
/* -------------------------------------------------------------------------- */
template<typename T, template<class> class ret_type>
class DumperIOHelper::material_strain_field_iterator :
public generic_internal_material_field_iterator<T, ret_type, StrainPaddingHelper,
material_strain_field_iterator> {
public:
typedef generic_internal_material_field_iterator<T, ret_type, StrainPaddingHelper,
DumperIOHelper::material_strain_field_iterator> parent;
typedef typename parent::it_type it_type;
typedef typename parent::data_type data_type;
typedef typename parent::return_type return_type;
typedef typename parent::field_type field_type;
typedef typename parent::internal_iterator internal_iterator;
typedef typename Vector<T>::template const_iterator< ret_type<T> > internal_material_iterator;
public:
material_strain_field_iterator(const field_type & element_material,
UInt n,
const typename field_type::type_iterator & t_it,
const typename field_type::type_iterator & t_it_end,
const internal_iterator & it,
ElementType element_type,
const GhostType ghost_type = _not_ghost) :
parent(element_material, n, t_it, t_it_end,
it, element_type, ghost_type) { }
};
/* -------------------------------------------------------------------------- */
template<typename T, template<class> class ret_type,
template<typename, template<class> class> class iterator_type>
class DumperIOHelper::InternalMaterialField : public GenericQuadraturePointsField<UInt,
- iterator_type<T, ret_type>,
- ret_type> {
+ iterator_type<T, ret_type>,
+ ret_type> {
+
+ typedef GenericQuadraturePointsField<UInt, iterator_type<T, ret_type>, ret_type> parent;
public:
- typedef iterator_type<T, ret_type> iterator;
-private:
- typedef GenericQuadraturePointsField<UInt, iterator, ret_type> parent;
+ typedef typename parent::iterator iterator;
public:
/* ------------------------------------------------------------------------ */
InternalMaterialField(const SolidMechanicsModel & model,
const std::string & field_id,
UInt spatial_dimension = 0,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) :
parent(model.getFEM(), model.getElementIndexByMaterial(), 0, spatial_dimension, ghost_type, element_kind),
model(model), field_id(field_id) {
init();
}
InternalMaterialField(const SolidMechanicsModel & model,
const std::string & field_id,
UInt n,
UInt spatial_dimension = 0,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) :
parent(model.getFEM(), model.getElementIndexByMaterial(), n, spatial_dimension, ghost_type, element_kind),
model(model), field_id(field_id) {
init();
}
virtual iterator begin() {
iterator it = parent::begin();
it.setItSize(2,1);
it.setModel(model);
it.setFieldID(field_id);
return it;
}
virtual iterator end () {
iterator it = parent::end();
it.setItSize(2,1);
it.setModel(model);
it.setFieldID(field_id);
return it;
}
protected:
void init() {
typename ByElementTypeVector<UInt>::type_iterator tit =
this->field.firstType(this->spatial_dimension, this->ghost_type, this->element_kind);
typename ByElementTypeVector<UInt>::type_iterator end =
this->field.lastType (this->spatial_dimension, this->ghost_type, this->element_kind);
UInt nb_materials = model.getNbMaterials();
bool found = false;
for(;tit != end; ++tit) {
// UInt nb_quad = this->fem.getNbQuadraturePoints(*tit);
for (UInt ma = 0; ma < nb_materials; ++ma) {
const Material & mat = model.getMaterial(ma);
try {
const Vector<Real> & vect __attribute__((unused)) = mat.getVector(field_id, *tit, this->ghost_type);
found = true;
} catch (...) { };
}
}
if(!found) AKANTU_DEBUG_ERROR("The field " << field_id << " does not exists in any materials");
}
virtual UInt getNbDataPerElem(__attribute__((unused)) const ElementType & type) { return 1; }
protected:
const SolidMechanicsModel & model;
std::string field_id;
};
diff --git a/src/io/dumper/dumper_iohelper_tmpl_quadrature_points_field.hh b/src/io/dumper/dumper_iohelper_tmpl_quadrature_points_field.hh
index ff7279aa8..5c54e1ea7 100644
--- a/src/io/dumper/dumper_iohelper_tmpl_quadrature_points_field.hh
+++ b/src/io/dumper/dumper_iohelper_tmpl_quadrature_points_field.hh
@@ -1,117 +1,172 @@
/**
* @file dumper_iohelper_tmpl_quadrature_points_field.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Fri Oct 26 21:52:40 2012
*
* @brief Description of quadrature points fields
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
__END_AKANTU__
#include "fem.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<typename i_type, typename d_type,
template<typename> class ret_type, class daughter>
-class DumperIOHelper::quadrature_point_iterator : public element_iterator<i_type, d_type,
- ret_type, daughter> {
+class DumperIOHelper::generic_quadrature_point_iterator : public element_iterator<i_type, d_type,
+ ret_type, daughter> {
public:
typedef element_iterator<i_type, d_type, ret_type, daughter> parent;
typedef typename parent::it_type it_type;
typedef typename parent::data_type data_type;
typedef typename parent::return_type return_type;
typedef typename parent::field_type field_type;
typedef typename parent::internal_iterator internal_iterator;
public:
- quadrature_point_iterator(const field_type & field,
- UInt n,
- const typename field_type::type_iterator & t_it,
- const typename field_type::type_iterator & t_it_end,
- const internal_iterator & it,
- ElementType element_type,
- const GhostType ghost_type = _not_ghost) :
+ generic_quadrature_point_iterator(const field_type & field,
+ UInt n,
+ const typename field_type::type_iterator & t_it,
+ const typename field_type::type_iterator & t_it_end,
+ const internal_iterator & it,
+ ElementType element_type,
+ const GhostType ghost_type = _not_ghost) :
parent(field, n, t_it, t_it_end, it, element_type, ghost_type) { }
void setFEM(const FEM & fem) { this->fem = &fem; }
protected:
virtual UInt getNbDataPerElem(const ElementType & type) { return fem->getNbQuadraturePoints(type); }
protected:
const FEM * fem;
};
+/* -------------------------------------------------------------------------- */
+template<typename T, template<typename> class ret_type>
+class DumperIOHelper::quadrature_point_iterator : public generic_quadrature_point_iterator<T, T,
+ ret_type, quadrature_point_iterator<T, ret_type> > {
+ public:
+ typedef generic_quadrature_point_iterator<T, T, ret_type, quadrature_point_iterator> parent;
+ typedef typename parent::it_type it_type;
+ typedef typename parent::data_type data_type;
+ typedef typename parent::return_type return_type;
+ typedef typename parent::field_type field_type;
+ typedef typename parent::internal_iterator internal_iterator;
+public:
+ quadrature_point_iterator(const field_type & field,
+ UInt n,
+ const typename field_type::type_iterator & t_it,
+ const typename field_type::type_iterator & t_it_end,
+ const internal_iterator & it,
+ ElementType element_type,
+ const GhostType ghost_type = _not_ghost) :
+ parent(field, n, t_it, t_it_end, it, element_type, ghost_type) { }
+
+ return_type operator*() {
+ UInt nb_data = this->fem->getNbQuadraturePoints(*this->tit);
+ return PaddingHelper<T, ret_type>::pad(*this->vit, this->padding_m, this->padding_n, nb_data);
+ }
+};
/* -------------------------------------------------------------------------- */
/* Fields type description */
/* -------------------------------------------------------------------------- */
-template<typename T, class iterator_type, template<class> class ret_type>
+template<typename T,
+ class iterator_type,
+ template<class> class ret_type>
class DumperIOHelper::GenericQuadraturePointsField : public GenericElementalField<T, iterator_type, ret_type> {
- typedef GenericElementalField<T, iterator_type, ret_type> parent;
public:
+ typedef iterator_type iterator;
+ typedef GenericElementalField<T, iterator, ret_type> parent;
+
GenericQuadraturePointsField(const FEM & fem,
const ByElementTypeVector<T> & field,
UInt spatial_dimension = 0,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) :
parent(field, spatial_dimension,
ghost_type, element_kind), fem(fem) { }
GenericQuadraturePointsField(const FEM & fem,
const ByElementTypeVector<T> & field,
UInt n,
UInt spatial_dimension = 0,
GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) :
parent(field, n, spatial_dimension,
ghost_type, element_kind), fem(fem) { }
- typedef iterator_type iterator;
-
/* ------------------------------------------------------------------------ */
virtual iterator begin() {
iterator it = parent::begin();
it.setFEM(fem);
return it;
}
virtual iterator end () {
iterator it = parent::end();
it.setFEM(fem);
return it;
}
const FEM & getFEM() { return fem; }
protected:
virtual UInt getNbDataPerElem(const ElementType & type) { return fem.getNbQuadraturePoints(type); }
protected:
const FEM & fem;
};
+
+/* -------------------------------------------------------------------------- */
+template<typename T,
+ template<class> class ret_type,
+ template<typename, template<class> class> class iterator_type>
+class DumperIOHelper::QuadraturePointsField : public GenericQuadraturePointsField<T, iterator_type<T, ret_type>, ret_type> {
+public:
+ typedef iterator_type<T, ret_type> iterator;
+ typedef GenericQuadraturePointsField<T, iterator, ret_type> parent;
+
+ QuadraturePointsField(const FEM & fem,
+ const ByElementTypeVector<T> & field,
+ UInt spatial_dimension = 0,
+ GhostType ghost_type = _not_ghost,
+ ElementKind element_kind = _ek_not_defined) :
+ parent(fem, field, spatial_dimension,
+ ghost_type, element_kind) { }
+
+ QuadraturePointsField(const FEM & fem,
+ const ByElementTypeVector<T> & field,
+ UInt n,
+ UInt spatial_dimension = 0,
+ GhostType ghost_type = _not_ghost,
+ ElementKind element_kind = _ek_not_defined) :
+ parent(fem, field, n, spatial_dimension,
+ ghost_type, element_kind) { }
+};
diff --git a/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc
index ef3878500..e526a1ec4 100644
--- a/src/model/heat_transfer/heat_transfer_model.cc
+++ b/src/model/heat_transfer/heat_transfer_model.cc
@@ -1,782 +1,813 @@
/**
* @file heat_transfer_model.cc
*
* @author Rui Wang <rui.wang@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
*
* @date Sun May 01 19:14:43 2011
*
* @brief Implementation of HeatTransferModel class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "heat_transfer_model.hh"
#include "aka_math.hh"
#include "aka_common.hh"
#include "fem_template.hh"
#include "mesh.hh"
#include "static_communicator.hh"
#include "parser.hh"
#include "generalized_trapezoidal.hh"
// #include "sparse_matrix.hh"
// #include "solver.hh"
// #ifdef AKANTU_USE_MUMPS
// #include "solver_mumps.hh"
// #endif
/* -------------------------------------------------------------------------- */
-
__BEGIN_AKANTU__
+#ifdef AKANTU_USE_IOHELPER
+# include "dumper_iohelper_tmpl_homogenizing_field.hh"
+#endif
+
/* -------------------------------------------------------------------------- */
HeatTransferModel::HeatTransferModel(Mesh & mesh,
- UInt dim,
- const ID & id,
- const MemoryID & memory_id) :
- Model(mesh, id, memory_id),
+ UInt dim,
+ const ID & id,
+ const MemoryID & memory_id) :
+ Model(mesh, id, memory_id), Dumpable<DumperParaview>(id),
integrator(new ForwardEuler()),
- spatial_dimension(dim),
+ spatial_dimension(dim != 0 ? dim : mesh.getSpatialDimension()),
temperature_gradient ("temperature_gradient", id),
temperature_on_qpoints ("temperature_on_qpoints", id),
conductivity_on_qpoints ("conductivity_on_qpoints", id),
k_gradt_on_qpoints ("k_gradt_on_qpoints", id),
int_bt_k_gT ("int_bt_k_gT", id),
bt_k_gT ("bt_k_gT", id),
+ conductivity(spatial_dimension, spatial_dimension),
thermal_energy ("thermal_energy", id)
{
AKANTU_DEBUG_IN();
createSynchronizerRegistry(this);
- if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension();
-
std::stringstream sstr; sstr << id << ":fem";
registerFEMObject<MyFEMType>(sstr.str(), mesh,spatial_dimension);
this->temperature= NULL;
-
this->residual = NULL;
this->boundary = NULL;
+
this->conductivity_variation = 0.0;
- this->t_ref = 0;
+ this->T_ref = 0;
+
+ addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::initModel() {
getFEM().initShapeFunctions(_not_ghost);
getFEM().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::initParallel(MeshPartition * partition,
- DataAccessor * data_accessor) {
+ DataAccessor * data_accessor) {
AKANTU_DEBUG_IN();
if (data_accessor == NULL) data_accessor = this;
Synchronizer & synch_parallel = createParallelSynch(partition,data_accessor);
synch_registry->registerSynchronizer(synch_parallel, _gst_htm_capacity);
synch_registry->registerSynchronizer(synch_parallel, _gst_htm_temperature);
synch_registry->registerSynchronizer(synch_parallel, _gst_htm_gradient_temperature);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::initPBC() {
AKANTU_DEBUG_IN();
-
+
Model::initPBC();
PBCSynchronizer * synch = new PBCSynchronizer(pbc_pair);
-
+
synch_registry->registerSynchronizer(*synch, _gst_htm_capacity);
synch_registry->registerSynchronizer(*synch, _gst_htm_temperature);
changeLocalEquationNumberforPBC(pbc_pair,1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::initVectors() {
AKANTU_DEBUG_IN();
UInt nb_nodes = getFEM().getMesh().getNbNodes();
std::stringstream sstr_temp; sstr_temp << id << ":temperature";
std::stringstream sstr_temp_rate; sstr_temp_rate << id << ":temperature_rate";
std::stringstream sstr_inc; sstr_inc << id << ":increment";
std::stringstream sstr_ext_flx; sstr_ext_flx << id << ":external_flux";
std::stringstream sstr_residual; sstr_residual << id << ":residual";
std::stringstream sstr_lump; sstr_lump << id << ":lumped";
std::stringstream sstr_boun; sstr_boun << id << ":boundary";
temperature = &(alloc<Real>(sstr_temp.str(), nb_nodes, 1, REAL_INIT_VALUE));
temperature_rate = &(alloc<Real>(sstr_temp_rate.str(), nb_nodes, 1, REAL_INIT_VALUE));
increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, 1, REAL_INIT_VALUE));
external_heat_rate = &(alloc<Real>(sstr_ext_flx.str(), nb_nodes, 1, REAL_INIT_VALUE));
residual = &(alloc<Real>(sstr_residual.str(), nb_nodes, 1, REAL_INIT_VALUE));
capacity_lumped = &(alloc<Real>(sstr_lump.str(), nb_nodes, 1, REAL_INIT_VALUE));
boundary = &(alloc<bool>(sstr_boun.str(), nb_nodes, 1, false));
Mesh::ConnectivityTypeList::const_iterator it;
/* -------------------------------------------------------------------------- */
// byelementtype vectors
-
getFEM().getMesh().initByElementTypeVector(temperature_on_qpoints,
- 1,
- spatial_dimension);
+ 1,
+ spatial_dimension);
getFEM().getMesh().initByElementTypeVector(temperature_gradient,
- spatial_dimension,
- spatial_dimension);
+ spatial_dimension,
+ spatial_dimension);
getFEM().getMesh().initByElementTypeVector(conductivity_on_qpoints,
- spatial_dimension*spatial_dimension,
- spatial_dimension);
+ spatial_dimension*spatial_dimension,
+ spatial_dimension);
getFEM().getMesh().initByElementTypeVector(k_gradt_on_qpoints,
- spatial_dimension,
- spatial_dimension);
+ spatial_dimension,
+ spatial_dimension);
getFEM().getMesh().initByElementTypeVector(bt_k_gT,
- 1,
- spatial_dimension,true);
+ 1,
+ spatial_dimension,true);
getFEM().getMesh().initByElementTypeVector(int_bt_k_gT,
- 1,
- spatial_dimension,true);
-
+ 1,
+ spatial_dimension,true);
+
getFEM().getMesh().initByElementTypeVector(thermal_energy,
- 1,
- spatial_dimension);
+ 1,
+ spatial_dimension);
for(UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType) g;
const Mesh::ConnectivityTypeList & type_list =
getFEM().getMesh().getConnectivityTypeList(gt);
for(it = type_list.begin(); it != type_list.end(); ++it) {
if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
UInt nb_element = getFEM().getMesh().getNbElement(*it, gt);
UInt nb_quad_points = this->getFEM().getNbQuadraturePoints(*it, gt) * nb_element;
temperature_on_qpoints(*it, gt).resize(nb_quad_points);
temperature_on_qpoints(*it, gt).clear();
temperature_gradient(*it, gt).resize(nb_quad_points);
temperature_gradient(*it, gt).clear();
conductivity_on_qpoints(*it, gt).resize(nb_quad_points);
conductivity_on_qpoints(*it, gt).clear();
k_gradt_on_qpoints(*it, gt).resize(nb_quad_points);
k_gradt_on_qpoints(*it, gt).clear();
bt_k_gT(*it, gt).resize(nb_quad_points);
bt_k_gT(*it, gt).clear();
int_bt_k_gT(*it, gt).resize(nb_element);
int_bt_k_gT(*it, gt).clear();
thermal_energy(*it, gt).resize(nb_element);
thermal_energy(*it, gt).clear();
}
}
/* -------------------------------------------------------------------------- */
dof_synchronizer = new DOFSynchronizer(getFEM().getMesh(),1);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
HeatTransferModel::~HeatTransferModel()
{
AKANTU_DEBUG_IN();
- delete [] conductivity;
if(dof_synchronizer) delete dof_synchronizer;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::assembleCapacityLumped(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
FEM & fem = getFEM();
const Mesh::ConnectivityTypeList & type_list
= fem.getMesh().getConnectivityTypeList(ghost_type);
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it)
{
if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
UInt nb_element = getFEM().getMesh().getNbElement(*it,ghost_type);
UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, ghost_type);
Vector<Real> rho_1 (nb_element * nb_quadrature_points,1, capacity * density);
fem.assembleFieldLumped(rho_1,1,*capacity_lumped,
- dof_synchronizer->getLocalDOFEquationNumbers(),
- *it, ghost_type);
+ dof_synchronizer->getLocalDOFEquationNumbers(),
+ *it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::assembleCapacityLumped() {
AKANTU_DEBUG_IN();
capacity_lumped->clear();
assembleCapacityLumped(_not_ghost);
assembleCapacityLumped(_ghost);
getSynchronizerRegistry().synchronize(akantu::_gst_htm_capacity);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::updateResidual() {
AKANTU_DEBUG_IN();
/// @f$ r = q_{ext} - q_{int} - C \dot T @f$
-
// start synchronization
synch_registry->asynchronousSynchronize(_gst_htm_temperature);
// finalize communications
synch_registry->waitEndSynchronize(_gst_htm_temperature);
//clear the array
/// first @f$ r = q_{ext} @f$
// residual->clear();
residual->copy(*external_heat_rate);
/// then @f$ r -= q_{int} @f$
// update the not ghost ones
updateResidual(_not_ghost);
// update for the received ghosts
updateResidual(_ghost);
/// finally @f$ r -= C \dot T @f$
// lumped C
UInt nb_nodes = temperature_rate->getSize();
UInt nb_degree_of_freedom = temperature_rate->getNbComponent();
Real * capacity_val = capacity_lumped->values;
Real * temp_rate_val = temperature_rate->values;
Real * res_val = residual->values;
bool * boundary_val = boundary->values;
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if(!(*boundary_val)) {
*res_val -= *capacity_val * *temp_rate_val;
}
boundary_val++;
res_val++;
capacity_val++;
temp_rate_val++;
}
#ifndef AKANTU_NDEBUG
getSynchronizerRegistry().synchronize(akantu::_gst_htm_gradient_temperature);
#endif
solveExplicitLumped();
AKANTU_DEBUG_OUT();
}
-
/* -------------------------------------------------------------------------- */
void HeatTransferModel::computeConductivityOnQuadPoints(const GhostType & ghost_type) {
const Mesh::ConnectivityTypeList & type_list =
this->getFEM().getMesh().getConnectivityTypeList(ghost_type);
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
+ Vector<Real> & temperature_interpolated = temperature_on_qpoints(*it, ghost_type);
+
//compute the temperature on quadrature points
this->getFEM().interpolateOnQuadraturePoints(*temperature,
- temperature_on_qpoints(*it, ghost_type),
- 1 ,*it,ghost_type);
+ temperature_interpolated,
+ 1 ,*it,ghost_type);
- Vector<Real>::iterator<types::RMatrix> c_iterator =
- conductivity_on_qpoints(*it,ghost_type).begin(spatial_dimension,spatial_dimension);
+ Vector<Real>::iterator<types::RMatrix> C_it =
+ conductivity_on_qpoints(*it, ghost_type).begin(spatial_dimension, spatial_dimension);
+ Vector<Real>::iterator<types::RMatrix> C_end =
+ conductivity_on_qpoints(*it, ghost_type).end(spatial_dimension, spatial_dimension);
- Real * T_val = temperature_on_qpoints(*it, ghost_type).values;
+ Vector<Real>::iterator<Real> T_it = temperature_interpolated.begin();
- UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, ghost_type);
- UInt nb_element = getFEM().getMesh().getNbElement(*it,ghost_type);
- for (UInt el = 0; el < nb_element; ++el) {
- for (UInt q = 0; q < nb_quadrature_points; ++q) {
- Real * c_onquad_val = c_iterator->storage();
- for (UInt i = 0; i < spatial_dimension*spatial_dimension; ++i) {
- // if (i == 0 && *T_val - t_ref > 0)
- // AKANTU_DEBUG_INFO("conductivity is " << conductivity[i]
- // << " variation is " << conductivity_variation
- // << " temp increase is " << *T_val - t_ref
- // << " new conductivity is " << conductivity[i] + conductivity_variation* (*T_val - t_ref));
- c_onquad_val[i] = conductivity[i] + conductivity_variation* (*T_val - t_ref);
- }
- ++T_val;
- ++c_iterator;
- }
+ for (;C_it != C_end; ++C_it, ++T_it) {
+ types::Matrix<Real> & C = *C_it;
+ Real & T = *T_it;
+ C = conductivity;
+
+ types::Matrix<Real> variation(spatial_dimension, spatial_dimension, conductivity_variation * (T - T_ref));
+ C += conductivity_variation;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::computeKgradT(const GhostType & ghost_type) {
computeConductivityOnQuadPoints(ghost_type);
const Mesh::ConnectivityTypeList & type_list =
this->getFEM().getMesh().getConnectivityTypeList(ghost_type);
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
+ const ElementType & type = *it;
if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
- UInt nb_element = getFEM().getMesh().getNbElement(*it,ghost_type);
- UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, ghost_type);
-
-
Vector<Real> & gradient = temperature_gradient(*it, ghost_type);
- gradient.resize(nb_element * nb_quadrature_points);
-
this->getFEM().gradientOnQuadraturePoints(*temperature,
- gradient,
- 1 ,*it,ghost_type);
+ gradient,
+ 1 ,*it,ghost_type);
- Vector<Real>::iterator<types::RMatrix> c_iterator =
- conductivity_on_qpoints(*it,ghost_type).begin(spatial_dimension,spatial_dimension);
+ Vector<Real>::iterator<types::RMatrix> C_it =
+ conductivity_on_qpoints(*it, ghost_type).begin(spatial_dimension, spatial_dimension);
- Vector<Real>::iterator<types::RVector> gT_iterator =
- temperature_gradient(*it,ghost_type).begin(spatial_dimension);
+ Vector<Real>::iterator<types::RVector> BT_it = gradient.begin(spatial_dimension);
- Vector<Real>::iterator<types::RVector> k_gT_iterator =
- k_gradt_on_qpoints(*it,ghost_type).begin(spatial_dimension);
+ Vector<Real>::iterator<types::RVector> k_BT_it =
+ k_gradt_on_qpoints(type, ghost_type).begin(spatial_dimension);
+ Vector<Real>::iterator<types::RVector> k_BT_end =
+ k_gradt_on_qpoints(type, ghost_type).end(spatial_dimension);
- for (UInt el = 0; el < nb_element; ++el) {
- for (UInt i = 0; i < nb_quadrature_points; ++i){
- Math::matrixt_vector(spatial_dimension, spatial_dimension,
- c_iterator->storage(),
- gT_iterator->storage(),
- k_gT_iterator->storage());
- ++c_iterator;
- ++gT_iterator;
- ++k_gT_iterator;
- }
+ for (;k_BT_it != k_BT_end; ++k_BT_it, ++BT_it, ++C_it) {
+ types::Vector<Real> & k_BT = *k_BT_it;
+ types::Vector<Real> & BT = *BT_it;
+ types::Matrix<Real> & C = *C_it;
+
+ k_BT.mul<false>(C, BT);
}
}
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::updateResidual(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
const Mesh::ConnectivityTypeList & type_list =
this->getFEM().getMesh().getConnectivityTypeList(ghost_type);
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
Vector<Real> & shapes_derivatives =
const_cast<Vector<Real> &>(getFEM().getShapesDerivatives(*it,ghost_type));
- UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, ghost_type);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
- UInt nb_element = getFEM().getMesh().getNbElement(*it,ghost_type);
// compute k \grad T
computeKgradT(ghost_type);
- Vector<Real>::iterator<types::RVector> k_gT_iterator =
+ Vector<Real>::iterator<types::RVector> k_BT_it =
k_gradt_on_qpoints(*it,ghost_type).begin(spatial_dimension);
- Vector<Real>::iterator<types::RMatrix> shapesd_iterator =
- shapes_derivatives.begin(nb_nodes_per_element,spatial_dimension);
+ Vector<Real>::iterator<types::RMatrix> B_it =
+ shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element);
- Vector<Real>::iterator<types::RVector> bt_k_gT_iterator =
+ Vector<Real>::iterator<types::RVector> Bt_k_BT_it =
bt_k_gT(*it,ghost_type).begin(nb_nodes_per_element);
+ Vector<Real>::iterator<types::RVector> Bt_k_BT_end =
+ bt_k_gT(*it,ghost_type).end(nb_nodes_per_element);
- // UInt bt_k_gT_size = nb_nodes_per_element;
- // Vector<Real> * bt_k_gT =
- // new Vector <Real> (nb_quadrature_points*nb_element, bt_k_gT_size);
- // Real * bt_k_gT_val = bt_k_gT->values;
-
- for (UInt el = 0; el < nb_element; ++el) {
- for (UInt i = 0; i < nb_quadrature_points; ++i) {
- Math::matrix_vector(nb_nodes_per_element, spatial_dimension,
- shapesd_iterator->storage(),
- k_gT_iterator->storage(),
- bt_k_gT_iterator->storage());
+ for (;Bt_k_BT_it != Bt_k_BT_end; ++Bt_k_BT_it, ++B_it, ++k_BT_it) {
+ types::Vector<Real> & k_BT = *k_BT_it;
+ types::Vector<Real> & Bt_k_BT = *Bt_k_BT_it;
+ types::Matrix<Real> & B = *B_it;
- ++shapesd_iterator;
- ++k_gT_iterator;
- ++bt_k_gT_iterator;
- }
+ Bt_k_BT.mul<true>(B, k_BT);
}
this->getFEM().integrate(bt_k_gT(*it,ghost_type),
- int_bt_k_gT(*it,ghost_type),
- nb_nodes_per_element, *it,ghost_type);
+ int_bt_k_gT(*it,ghost_type),
+ nb_nodes_per_element, *it,ghost_type);
this->getFEM().assembleVector(int_bt_k_gT(*it,ghost_type), *residual,
- dof_synchronizer->getLocalDOFEquationNumbers(),
- 1, *it,ghost_type,NULL,-1);
- //delete q_e;
-
+ dof_synchronizer->getLocalDOFEquationNumbers(),
+ 1, *it,ghost_type,NULL,-1);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::solveExplicitLumped() {
AKANTU_DEBUG_IN();
UInt nb_nodes = increment->getSize();
UInt nb_degree_of_freedom = increment->getNbComponent();
Real * capa_val = capacity_lumped->values;
Real * res_val = residual->values;
bool * boundary_val = boundary->values;
Real * inc = increment->values;
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if(!(*boundary_val)) {
*inc = (*res_val / *capa_val);
}
res_val++;
boundary_val++;
inc++;
capa_val++;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::explicitPred() {
AKANTU_DEBUG_IN();
integrator->integrationSchemePred(time_step,
- *temperature,
- *temperature_rate,
- *boundary);
+ *temperature,
+ *temperature_rate,
+ *boundary);
UInt nb_nodes = temperature->getSize();
UInt nb_degree_of_freedom = temperature->getNbComponent();
Real * temp = temperature->values;
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n, ++temp)
if(*temp < 0.) *temp = 0.;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::explicitCorr() {
AKANTU_DEBUG_IN();
integrator->integrationSchemeCorrTempRate(time_step,
- *temperature,
- *temperature_rate,
- *boundary,
- *increment);
+ *temperature,
+ *temperature_rate,
+ *boundary,
+ *increment);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real HeatTransferModel::getStableTimeStep()
{
AKANTU_DEBUG_IN();
- Real conductivitymax = -std::numeric_limits<Real>::max();
Real el_size;
Real min_el_size = std::numeric_limits<Real>::max();
+ Real conductivitymax = conductivity(0, 0);
+
+ //get the biggest parameter from k11 until k33//
+ for(UInt i = 0; i < spatial_dimension; i++)
+ for(UInt j = 0; j < spatial_dimension; j++)
+ conductivitymax = std::max(conductivity(i, j), conductivitymax);
+
const Mesh::ConnectivityTypeList & type_list =
getFEM().getMesh().getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
if(getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension)
continue;
UInt nb_nodes_per_element = getFEM().getMesh().getNbNodesPerElement(*it);
Vector<Real> coord(0, nb_nodes_per_element*spatial_dimension);
FEM::extractNodalToElementField(getFEM().getMesh(), getFEM().getMesh().getNodes(),
coord, *it, _not_ghost);
Vector<Real>::iterator< types::Matrix<Real> > el_coord = coord.begin(spatial_dimension, nb_nodes_per_element);
UInt nb_element = getFEM().getMesh().getNbElement(*it);
for (UInt el = 0; el < nb_element; ++el, ++el_coord) {
- el_size = getFEM().getElementInradius(*el_coord, *it);
- min_el_size = std::min(min_el_size, el_size);
+ el_size = getFEM().getElementInradius(*el_coord, *it);
+ min_el_size = std::min(min_el_size, el_size);
}
AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size
- << " and the max conductivity is : "
- << conductivitymax);
- }
-
- //get the biggest parameter from k11 until k33//
- for(UInt i = 0; i < spatial_dimension * spatial_dimension; i++) {
- if(conductivity[i] > conductivitymax)
- conductivitymax=conductivity[i];
+ << " and the max conductivity is : "
+ << conductivitymax);
}
Real min_dt = 2 * min_el_size * min_el_size * density
- * capacity/conductivitymax;
+ * capacity / conductivitymax;
StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::readMaterials(const std::string & filename) {
Parser parser;
parser.open(filename);
std::string opt_param;
std::string mat_type = parser.getNextSection("heat", opt_param);
- if (mat_type != ""){
+ if (mat_type != "") {
parser.readSection<HeatTransferModel>(*this);
}
else
AKANTU_DEBUG_ERROR("did not find any section with material info "
- <<"for heat conduction");
+ <<"for heat conduction");
}
/* -------------------------------------------------------------------------- */
bool HeatTransferModel::setParam(const std::string & key,
- const std::string & value) {
+ const std::string & value) {
std::stringstream str(value);
if (key == "conductivity") {
- conductivity = new Real[spatial_dimension * spatial_dimension];
for(UInt i = 0; i < 3; i++)
for(UInt j = 0; j < 3; j++) {
- if (i< spatial_dimension && j < spatial_dimension){
- str >> conductivity[i*spatial_dimension+j];
- AKANTU_DEBUG_INFO("Conductivity(" << i << "," << j << ") = "
- << conductivity[i*spatial_dimension+j]);
- } else {
- Real tmp;
- str >> tmp;
- }
+ if (i< spatial_dimension && j < spatial_dimension){
+ str >> conductivity(i, j);
+ AKANTU_DEBUG_INFO("Conductivity(" << i << "," << j << ") = "
+ << conductivity[i*spatial_dimension+j]);
+ } else {
+ Real tmp;
+ str >> tmp;
+ }
}
}
else if (key == "conductivity_variation") {
str >> conductivity_variation;
}
else if (key == "temperature_reference") {
- str >> t_ref;
+ str >> T_ref;
}
else if (key == "capacity"){
str >> capacity;
AKANTU_DEBUG_INFO("The capacity of the material is:" << capacity);
}
else if (key == "density"){
str >> density;
AKANTU_DEBUG_INFO("The density of the material is:" << density);
} else {
return false;
}
return true;
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::initFull(const std::string & material_file){
readMaterials(material_file);
//model initialization
initModel();
//initialize the vectors
initVectors();
temperature->clear();
temperature_rate->clear();
external_heat_rate->clear();
}
/* -------------------------------------------------------------------------- */
void HeatTransferModel::initFEMBoundary(bool create_surface) {
if(create_surface)
MeshUtils::buildFacets(getFEM().getMesh());
FEM & fem_boundary = getFEMBoundary();
fem_boundary.initShapeFunctions();
fem_boundary.computeNormalsOnControlPoints();
}
/* -------------------------------------------------------------------------- */
-
Real HeatTransferModel::computeThermalEnergyByNode() {
AKANTU_DEBUG_IN();
Real ethermal = 0.;
Vector<Real>::iterator<types::RVector> heat_rate_it =
residual->begin(residual->getNbComponent());
Vector<Real>::iterator<types::RVector> heat_rate_end =
residual->end(residual->getNbComponent());
UInt n = 0;
for(;heat_rate_it != heat_rate_end; ++heat_rate_it, ++n) {
Real heat = 0;
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
-
+
types::RVector & heat_rate = *heat_rate_it;
for (UInt i = 0; i < heat_rate.size(); ++i) {
if (count_node)
heat += heat_rate[i] * time_step;
}
ethermal += heat;
}
StaticCommunicator::getStaticCommunicator().allReduce(ðermal, 1, _so_sum);
AKANTU_DEBUG_OUT();
return ethermal;
}
/* -------------------------------------------------------------------------- */
+template<class iterator>
+void HeatTransferModel::getThermalEnergy(iterator Eth,
+ Vector<Real>::const_iterator<Real> T_it,
+ Vector<Real>::const_iterator<Real> T_end) const {
+ for(;T_it != T_end; ++T_it, ++Eth) {
+ *Eth = capacity * density * *T_it;
+ }
+}
-void HeatTransferModel::computeThermalEnergyByElement() {
+/* -------------------------------------------------------------------------- */
+Real HeatTransferModel::getThermalEnergy(const ElementType & type, UInt index) {
AKANTU_DEBUG_IN();
- Mesh & mesh = getFEM().getMesh();
- Mesh::type_iterator it = mesh.firstType(spatial_dimension);
- Mesh::type_iterator last_type = mesh.lastType(spatial_dimension);
- for(; it != last_type; ++it) {
- if(!thermal_energy.exists(*it, _not_ghost)) {
- UInt nb_element = getFEM().getMesh().getNbElement(*it, _not_ghost);
- UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, _not_ghost);
+ UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(type);
+ types::Vector<Real> Eth_on_quarature_points(nb_quadrature_points);
- thermal_energy.alloc(nb_element * nb_quadrature_points, 1,
- *it, _not_ghost);
- }
+ Vector<Real>::iterator<Real> T_it = this->temperature_on_qpoints(type).begin();
+ T_it += index * nb_quadrature_points;
- Real * ethermal = thermal_energy(*it, _not_ghost).storage();
- Vector<Real>::iterator<Real> temperature_it =
- this->temperature_on_qpoints(*it).begin();
- Vector<Real>::iterator<Real> temperature_end =
- this->temperature_on_qpoints(*it).end();
+ Vector<Real>::iterator<Real> T_end = T_it + nb_quadrature_points;
- for(;temperature_it != temperature_end; ++temperature_it, ++ ethermal) {
- *ethermal = capacity * density * *temperature_it;
- }
- }
- AKANTU_DEBUG_OUT();
-}
+ getThermalEnergy(Eth_on_quarature_points.storage(), T_it, T_end);
-Real HeatTransferModel::getEnergy(const std::string & id) {
- AKANTU_DEBUG_IN();
+ return getFEM().integrate(Eth_on_quarature_points, type, index);
+}
- //return computeThermalEnergyByNode();
+/* -------------------------------------------------------------------------- */
+Real HeatTransferModel::getThermalEnergy() {
+ Real Eth = 0;
- Real energy = 0.;
- computeThermalEnergyByElement();
- /// integrate the thermal energy for each type of element
Mesh & mesh = getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension);
for(; it != last_type; ++it) {
+ UInt nb_element = getFEM().getMesh().getNbElement(*it, _not_ghost);
+ UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(*it, _not_ghost);
+ Vector<Real> Eth_per_quad(nb_element * nb_quadrature_points, 1);
+
+ Vector<Real>::iterator<Real> T_it = this->temperature_on_qpoints(*it).begin();
+ Vector<Real>::iterator<Real> T_end = this->temperature_on_qpoints(*it).end();
+ getThermalEnergy(Eth_per_quad.begin(), T_it, T_end);
- energy += getFEM().integrate(thermal_energy(*it, _not_ghost), *it,
- _not_ghost);
+ Eth += getFEM().integrate(Eth, *it);
}
- /// reduction sum over all processors
+ return Eth;
+}
+
+/* -------------------------------------------------------------------------- */
+Real HeatTransferModel::getEnergy(const std::string & id) {
+ AKANTU_DEBUG_IN();
+ Real energy = 0;
+
+ if("thermal") energy = getThermalEnergy();
+
+ // reduction sum over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
-
Real HeatTransferModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index) {
AKANTU_DEBUG_IN();
- Real ethermal = 0.;
- types::RVector ethermal_on_quad_points(getFEM().getNbQuadraturePoints(type));
+ Real energy = 0.;
- Vector<Real>::iterator<Real> temperature_it =
- this->temperature_on_qpoints(type).begin();
- Vector<Real>::iterator<Real> temperature_end =
- this->temperature_on_qpoints(type).begin();
+ if("thermal") energy = getThermalEnergy(type, index);
- UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(type);
+ AKANTU_DEBUG_OUT();
+ return energy;
+}
- temperature_it += index*nb_quadrature_points;
- temperature_end += (index+1)*nb_quadrature_points;
+/* -------------------------------------------------------------------------- */
+#ifdef AKANTU_USE_IOHELPER
+#define IF_ADD_NODAL_FIELD(field, type, pad) \
+ if(field_id == BOOST_PP_STRINGIZE(field)) { \
+ DumperIOHelper::Field * f = \
+ new DumperIOHelper::NodalField<type>(*field); \
+ if(pad) f->setPadding(3); \
+ addDumpFieldToDumper(BOOST_PP_STRINGIZE(field), f); \
+ } else
+
+#define IF_ADD_ELEM_FIELD(field, type, pad) \
+ if(field_id == BOOST_PP_STRINGIZE(field)) { \
+ typedef DumperIOHelper::HomogenizedField<type, \
+ DumperIOHelper::QuadraturePointsField, \
+ DumperIOHelper::quadrature_point_iterator> Field; \
+ DumperIOHelper::Field * f = \
+ new Field(getFEM(), \
+ field, \
+ spatial_dimension, \
+ spatial_dimension, \
+ _not_ghost, \
+ _ek_regular); \
+ if(pad) f->setPadding(3); \
+ addDumpFieldToDumper(BOOST_PP_STRINGIZE(field), f); \
+ } else
+#endif
- Real * ethermal_quad = ethermal_on_quad_points.storage();
- for(;temperature_it != temperature_end; ++temperature_it, ++ ethermal_quad) {
- *ethermal_quad = capacity * density * *temperature_it;
+void HeatTransferModel::addDumpField(const std::string & field_id) {
+#ifdef AKANTU_USE_IOHELPER
+ IF_ADD_NODAL_FIELD(temperature , Real, false)
+ IF_ADD_NODAL_FIELD(temperature_rate , Real, false)
+ IF_ADD_NODAL_FIELD(external_heat_rate, Real, false)
+ IF_ADD_NODAL_FIELD(residual , Real, false)
+ IF_ADD_NODAL_FIELD(capacity_lumped , Real, false)
+ IF_ADD_NODAL_FIELD(boundary , bool, false)
+ IF_ADD_ELEM_FIELD (temperature_gradient, Real, false)
+ if(field_id == "partitions" ) {
+ addDumpFieldToDumper(field_id,
+ new DumperIOHelper::ElementPartitionField(mesh,
+ spatial_dimension,
+ _not_ghost,
+ _ek_regular));
+ } else {
+ AKANTU_DEBUG_ERROR("Field " << field_id << " does not exists in the model " << id);
}
-
- ethermal = getFEM().integrate(ethermal_on_quad_points, type, index);
+#endif
+}
- AKANTU_DEBUG_OUT();
- return ethermal;
+/* -------------------------------------------------------------------------- */
+void HeatTransferModel::addDumpFieldVector(const std::string & field_id) {
+#ifdef AKANTU_USE_IOHELPER
+ IF_ADD_ELEM_FIELD (temperature_gradient, Real, false)
+ {
+ AKANTU_DEBUG_ERROR("Field " << field_id << " does not exists in the model " << id
+ << " or is not dumpable as a vector");
+ }
+#endif
}
/* -------------------------------------------------------------------------- */
+void HeatTransferModel::addDumpFieldTensor(const std::string & field_id) {
+#ifdef AKANTU_USE_IOHELPER
+ AKANTU_DEBUG_ERROR("Field " << field_id << " does not exists in the model " << id
+ << " or is not dumpable as a Tensor");
+#endif
+}
__END_AKANTU__
diff --git a/src/model/heat_transfer/heat_transfer_model.hh b/src/model/heat_transfer/heat_transfer_model.hh
index 478cc8ac4..ecfa452b8 100644
--- a/src/model/heat_transfer/heat_transfer_model.hh
+++ b/src/model/heat_transfer/heat_transfer_model.hh
@@ -1,318 +1,336 @@
/**
* @file heat_transfer_model.hh
*
* @author Rui Wang <rui.wang@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Sun May 01 19:14:43 2011
*
* @brief Model of Heat Transfer
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_HEAT_TRANSFER_MODEL_HH__
#define __AKANTU_HEAT_TRANSFER_MODEL_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_memory.hh"
#include "model.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
+#include "dumpable.hh"
+
namespace akantu {
class IntegrationScheme1stOrder;
-// class Solver;
-// class SparseMatrix;
}
__BEGIN_AKANTU__
-class HeatTransferModel : public Model, public DataAccessor {
+class HeatTransferModel : public Model, public DataAccessor, public Dumpable<DumperParaview> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
typedef FEMTemplate<IntegratorGauss,ShapeLagrange> MyFEMType;
HeatTransferModel(UInt spatial_dimension,
const ID & id = "heat_transfer_model",
const MemoryID & memory_id = 0) ;
HeatTransferModel(Mesh & mesh,
UInt spatial_dimension = 0,
const ID & id = "heat_transfer_model",
const MemoryID & memory_id = 0);
virtual ~HeatTransferModel() ;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// generic function to initialize everything ready for explicit dynamics
void initFull(const std::string & material_file);
/// initialize the fem object of the boundary
void initFEMBoundary(bool create_surface = true);
/// set the parameters
bool setParam(const std::string & key, const std::string & value);
/// read one material file to instantiate all the materials
void readMaterials(const std::string & filename);
/// allocate all vectors
void initVectors();
/// register the tags associated with the parallel synchronizer
void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL);
/// initialize the model
void initModel();
/// init PBC synchronizer
void initPBC();
/// function to print the contain of the class
virtual void printself(__attribute__ ((unused)) std::ostream & stream,
__attribute__ ((unused)) int indent = 0) const {};
/* ------------------------------------------------------------------------ */
/* Methods for explicit */
/* ------------------------------------------------------------------------ */
public:
/// compute and get the stable time step
Real getStableTimeStep();
/// compute the heat flux
void updateResidual();
/// calculate the lumped capacity vector for heat transfer problem
void assembleCapacityLumped();
/// update the temperature from the temperature rate
void explicitPred();
/// update the temperature rate from the increment
void explicitCorr();
// /// initialize the heat flux
// void initializeResidual(Vector<Real> &temp);
// /// initialize temperature
// void initializeTemperature(Vector<Real> &temp);
private:
/// solve the system in temperature rate @f$C\delta \dot T = q_{n+1} - C \dot T_{n}@f$
void solveExplicitLumped();
/// compute the heat flux on ghost types
void updateResidual(const GhostType & ghost_type);
/// calculate the lumped capacity vector for heat transfer problem (w ghosttype)
void assembleCapacityLumped(const GhostType & ghost_type);
/// compute the conductivity tensor for each quadrature point in an array
void computeConductivityOnQuadPoints(const GhostType & ghost_type);
/// compute vector k \grad T for each quadrature point
void computeKgradT(const GhostType & ghost_type);
+ /// compute the thermal energy
+ Real computeThermalEnergyByNode();
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
inline UInt getNbDataForElements(const Vector<Element> & elements,
SynchronizationTag tag) const;
inline void packElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag) const;
inline void unpackElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag);
inline UInt getNbDataToPack(SynchronizationTag tag) const;
inline UInt getNbDataToUnpack(SynchronizationTag tag) const;
inline void packData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag) const;
inline void unpackData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag);
+ /* ------------------------------------------------------------------------ */
+ /* Dumpable interface */
+ /* ------------------------------------------------------------------------ */
+public:
+ virtual void addDumpField(const std::string & field_id);
+ virtual void addDumpFieldVector(const std::string & field_id);
+ virtual void addDumpFieldTensor(const std::string & field_id);
+
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
inline FEM & getFEMBoundary(std::string name = "");
/// get the dimension of the system space
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
/// get the current value of the time step
AKANTU_GET_MACRO(TimeStep, time_step, Real);
/// set the value of the time step
AKANTU_SET_MACRO(TimeStep, time_step, Real);
/// get the assembled heat flux
AKANTU_GET_MACRO(Residual, *residual, Vector<Real>&);
/// get the lumped capacity
AKANTU_GET_MACRO(CapacityLumped, * capacity_lumped, Vector<Real>&);
/// get the boundary vector
AKANTU_GET_MACRO(Boundary, * boundary, Vector<bool>&);
/// get the external heat rate vector
AKANTU_GET_MACRO(ExternalHeatRate, * external_heat_rate, Vector<Real>&);
/// get the temperature gradient
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureGradient, temperature_gradient, Real);
/// get the conductivity on q points
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConductivityOnQpoints, conductivity_on_qpoints, Real);
/// get the conductivity on q points
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureOnQpoints, temperature_on_qpoints, Real);
/// internal variables
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(KGradtOnQpoints, k_gradt_on_qpoints, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(IntBtKgT, int_bt_k_gT, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(BtKgT, bt_k_gT, Real);
/// get the temperature
AKANTU_GET_MACRO(Temperature, *temperature, Vector<Real> &);
/// get the temperature derivative
AKANTU_GET_MACRO(TemperatureRate, *temperature_rate, Vector<Real> &);
/// get the equation number Vector<Int>
AKANTU_GET_MACRO(EquationNumber, *equation_number, const Vector<Int> &);
- /// compute the thermal energy
- Real computeThermalEnergyByNode();
- /// get thermal energy by element
+ /// get the energy denominated by thermal
Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index);
- /// get thermal energy
+ /// get the energy denominated by thermal
Real getEnergy(const std::string & energy_id);
- /// compute thermal energy by element
- void computeThermalEnergyByElement();
+
+ /// get the thermal energy for a given element
+ Real getThermalEnergy(const ElementType & type, UInt index);
+ /// get the thermal energy for a given element
+ Real getThermalEnergy();
+
+protected:
+ /* ----------------------------------------------------------------------- */
+ template<class iterator>
+ void getThermalEnergy(iterator Eth,
+ Vector<Real>::const_iterator<Real> T_it,
+ Vector<Real>::const_iterator<Real> T_end) const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
IntegrationScheme1stOrder * integrator;
/// time step
Real time_step;
/// temperatures array
Vector<Real> * temperature;
/// temperatures derivatives array
Vector<Real> * temperature_rate;
/// increment array (@f$\delta \dot T@f$ or @f$\delta T@f$)
Vector<Real> * increment;
/// the spatial dimension
UInt spatial_dimension;
/// the density
Real density;
/// the speed of the changing temperature
ByElementTypeReal temperature_gradient;
/// temperature field on quadrature points
ByElementTypeReal temperature_on_qpoints;
/// conductivity tensor on quadrature points
ByElementTypeReal conductivity_on_qpoints;
/// vector k \grad T on quad points
ByElementTypeReal k_gradt_on_qpoints;
/// vector \int \grad N k \grad T
ByElementTypeReal int_bt_k_gT;
/// vector \grad N k \grad T
ByElementTypeReal bt_k_gT;
/// external flux vector
Vector<Real> * external_heat_rate;
/// residuals array
Vector<Real> * residual;
/// position of a dof in the K matrix
Vector<Int> * equation_number;
//lumped vector
Vector<Real> * capacity_lumped;
/// boundary vector
Vector<bool> * boundary;
//realtime
Real time;
///capacity
Real capacity;
//conductivity matrix
- Real* conductivity;
+ types::Matrix<Real> conductivity;
//linear variation of the conductivity (for temperature dependent conductivity)
Real conductivity_variation;
// reference temperature for the interpretation of temperature variation
- Real t_ref;
+ Real T_ref;
//the biggest parameter of conductivity matrix
Real conductivitymax;
/// thermal energy by element
ByElementTypeReal thermal_energy;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "heat_transfer_model_inline_impl.cc"
#endif
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const HeatTransferModel & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#endif /* __AKANTU_HEAT_TRANSFER_MODEL_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 289612d88..a3431433f 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,876 +1,876 @@
/**
* @file material_non_local_inline_impl.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Wed Aug 31 11:09:48 2011
*
* @brief Non-local inline implementation
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
__END_AKANTU__
/* -------------------------------------------------------------------------- */
#include "aka_types.hh"
#include "grid_synchronizer.hh"
#include "synchronizer_registry.hh"
#include "integrator.hh"
/* -------------------------------------------------------------------------- */
#include <fstream>
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt DIM, template <UInt> class WeightFunction>
MaterialNonLocal<DIM, WeightFunction>::MaterialNonLocal(SolidMechanicsModel & model,
const ID & id) :
Material(model, id), radius(100.), weight_func(NULL), spatial_grid(NULL),
update_weights(0), compute_stress_calls(0), is_creating_grid(false), grid_synchronizer(NULL) {
AKANTU_DEBUG_IN();
this->registerParam("radius" , radius , 100., _pat_readable , "Non local radius");
this->registerParam("UpdateWeights" , update_weights, 0U , _pat_modifiable, "Update weights frequency");
this->registerParam("Weight function", weight_func , _pat_internal);
this->is_non_local = true;
this->weight_func = new WeightFunction<DIM>(*this);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
MaterialNonLocal<spatial_dimension, WeightFunction>::~MaterialNonLocal() {
AKANTU_DEBUG_IN();
delete spatial_grid;
delete weight_func;
delete grid_synchronizer;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::initMaterial() {
AKANTU_DEBUG_IN();
// Material::initMaterial();
Mesh & mesh = this->model->getFEM().getMesh();
ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", id);
this->initInternalVector(quadrature_points_coordinates, spatial_dimension, true);
ByElementType<UInt> 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);
createCellList(quadrature_points_coordinates);
updatePairList(quadrature_points_coordinates);
#if not defined(AKANTU_NDEBUG)
neighbourhoodStatistics("material_non_local.stats");
#endif
- // cleanupExtraGhostElement(nb_ghost_protected);
+ cleanupExtraGhostElement(nb_ghost_protected);
weight_func->setRadius(radius);
weight_func->init();
computeWeights(quadrature_points_coordinates);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::cleanupExtraGhostElement(const ByElementType<UInt> & nb_ghost_protected) {
AKANTU_DEBUG_IN();
// Create list of element to keep
std::set<Element> relevant_ghost_element;
pair_type::const_iterator first_pair_types = existing_pairs[1].begin();
pair_type::const_iterator last_pair_types = existing_pairs[1].end();
for (; first_pair_types != last_pair_types; ++first_pair_types) {
ElementType type2 = first_pair_types->second;
GhostType ghost_type2 = _ghost;
UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2);
Vector<UInt> & elem_filter = element_filter(type2, ghost_type2);
const Vector<UInt> & pairs =
pair_list(first_pair_types->first, _not_ghost)(first_pair_types->second, ghost_type2);
Vector<UInt>::const_iterator< types::Vector<UInt> > first_pair = pairs.begin(2);
Vector<UInt>::const_iterator< types::Vector<UInt> > last_pair = pairs.end(2);
for(;first_pair != last_pair; ++first_pair) {
UInt _q2 = (*first_pair)(1);
QuadraturePoint q2(type2, elem_filter(_q2 / nb_quad2), _q2 % nb_quad2, ghost_type2);
relevant_ghost_element.insert(q2);
}
}
// Create list of element to remove and new numbering for element to keep
Mesh & mesh = this->model->getFEM().getMesh();
std::set<Element> 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;
element.ghost_type = _ghost;
for(; it != last_type; ++it) {
element.type = *it;
UInt nb_ghost_elem = 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);
Vector<UInt> & elem_filter = element_filter(*it, _ghost);
Vector<UInt> & new_numbering = remove_elem.getNewNumbering(*it, _ghost);
UInt ng = 0;
for (UInt g = 0; g < nb_ghost_elem; ++g) {
element.element = elem_filter(g);
if(element.element >= nb_ghost_elem_protected &&
(std::find(relevant_ghost_element.begin(),
relevant_ghost_element.end(),
element) == relevant_ghost_element.end())) {
ghost_to_erase.insert(element);
remove_elem.getList().push_back(element);
new_numbering(g) = UInt(-1);
} else {
new_numbering(g) = ng;
++ng;
}
}
}
mesh.sendEvent(remove_elem);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::createCellList(ByElementTypeReal & quadrature_points_coordinates) {
AKANTU_DEBUG_IN();
const Real safety_factor = 1.2; // for the cell grid spacing
Mesh & mesh = this->model->getFEM().getMesh();
mesh.computeBoundingBox();
Real lower_bounds[spatial_dimension];
Real upper_bounds[spatial_dimension];
mesh.getLocalLowerBounds(lower_bounds);
mesh.getLocalUpperBounds(upper_bounds);
types::Vector<Real> spacing(spatial_dimension, radius * safety_factor);
types::Vector<Real> center(spatial_dimension);
for (UInt i = 0; i < spatial_dimension; ++i) {
center(i) = (upper_bounds[i] + lower_bounds[i]) / 2.;
}
spatial_grid = new SpatialGrid<QuadraturePoint>(spatial_dimension, spacing, center);
computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost);
fillCellList(quadrature_points_coordinates, _not_ghost);
is_creating_grid = true;
SynchronizerRegistry & synch_registry = this->model->getSynchronizerRegistry();
std::stringstream sstr; sstr << id << ":grid_synchronizer";
grid_synchronizer = GridSynchronizer::createGridSynchronizer(mesh,
*spatial_grid,
sstr.str());
synch_registry.registerSynchronizer(*grid_synchronizer, _gst_mnl_for_average);
synch_registry.registerSynchronizer(*grid_synchronizer, _gst_mnl_weight);
is_creating_grid = false;
#if not defined(AKANTU_NDEBUG)
Mesh * mesh_tmp = new Mesh(spatial_dimension, "mnl_grid");
spatial_grid->saveAsMesh(*mesh_tmp);
std::stringstream sstr_grid;
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
Int prank = comm.whoAmI();
sstr_grid << "material_non_local_grid_" << std::setfill('0') << std::setw(4) << prank << ".msh";
mesh_tmp->write(sstr_grid.str());
delete mesh_tmp;
#endif
this->computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost);
fillCellList(quadrature_points_coordinates, _ghost);
#if not defined(AKANTU_NDEBUG)
mesh_tmp = new Mesh(spatial_dimension, "mnl_grid");
spatial_grid->saveAsMesh(*mesh_tmp);
sstr_grid.str(std::string());
sstr_grid << "material_non_local_grid_ghost_" << std::setfill('0') << std::setw(4) << prank << ".msh";
mesh_tmp->write(sstr_grid.str());
delete mesh_tmp;
#endif
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::fillCellList(const ByElementTypeReal & quadrature_points_coordinates,
const GhostType & ghost_type) {
Mesh & mesh = this->model->getFEM().getMesh();
QuadraturePoint q;
q.ghost_type = ghost_type;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType (spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
Vector<UInt> & elem_filter = element_filter(*it, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(*it, ghost_type);
const Vector<Real> & quads = quadrature_points_coordinates(*it, ghost_type);
q.type = *it;
Vector<Real>::const_iterator<types::RVector> 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.setPosition(*quad);
spatial_grid->insert(q, *quad);
++quad;
}
++elem;
}
}
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::updatePairList(const ByElementTypeReal & quadrature_points_coordinates) {
AKANTU_DEBUG_IN();
Mesh & mesh = this->model->getFEM().getMesh();
GhostType ghost_type = _not_ghost;
// generate the pair of neighbor depending of the cell_list
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
// Preparing datas
const Vector<Real> & quads = quadrature_points_coordinates(*it, ghost_type);
Vector<Real>::const_iterator<types::RVector> first_quad = quads.begin(spatial_dimension);
Vector<Real>::const_iterator<types::RVector> last_quad = quads.end(spatial_dimension);
ByElementTypeUInt & pairs = pair_list(ByElementTypeUInt("pairs", id, memory_id),
*it,
ghost_type);
ElementType current_element_type = _not_defined;
GhostType current_ghost_type = _casper;
UInt existing_pairs_num = 0;
Vector<UInt> * neighbors = NULL;
Vector<UInt> * element_index_material2 = NULL;
UInt my_num_quad = 0;
// loop over quad points
for(;first_quad != last_quad; ++first_quad, ++my_num_quad) {
SpatialGrid<QuadraturePoint>::CellID cell_id = spatial_grid->getCellID(*first_quad);
SpatialGrid<QuadraturePoint>::neighbor_cells_iterator first_neigh_cell =
spatial_grid->beginNeighborCells(cell_id);
SpatialGrid<QuadraturePoint>::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<QuadraturePoint>::Cell::iterator first_neigh_quad =
spatial_grid->beginCell(*first_neigh_cell);
SpatialGrid<QuadraturePoint>::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 & quad = *first_neigh_quad;
UInt nb_quad_per_elem =
this->model->getFEM().getNbQuadraturePoints(quad.type,
quad.ghost_type);
// little optimization to not search in the map at each quad points
if(quad.type != current_element_type ||
quad.ghost_type != current_ghost_type) {
current_element_type = quad.type;
current_ghost_type = quad.ghost_type;
existing_pairs_num = quad.ghost_type == _not_ghost ? 0 : 1;
if(!pairs.exists(current_element_type, current_ghost_type)) {
neighbors = &(pairs.alloc(0, 2,
current_element_type,
current_ghost_type));
} else {
neighbors = &(pairs(current_element_type,
current_ghost_type));
}
existing_pairs[existing_pairs_num].insert(std::pair<ElementType,
ElementType>(*it,
current_element_type));
element_index_material2 =
&(this->model->getElementIndexByMaterial(current_element_type,
current_ghost_type));
}
UInt neigh_num_quad =
(*element_index_material2)(quad.element) * nb_quad_per_elem +
quad.num_point;
const types::RVector & neigh_quad = quad.getPosition();
Real distance = first_quad->distance(neigh_quad);
if(distance <= radius &&
(current_ghost_type == _ghost ||
(current_ghost_type == _not_ghost && my_num_quad <= neigh_num_quad))) { // sotring only half lists
UInt pair[2];
pair[0] = my_num_quad;
pair[1] = neigh_num_quad;
neighbors->push_back(pair);
}
}
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::computeWeights(const ByElementTypeReal & quadrature_points_coordinates) {
AKANTU_DEBUG_IN();
GhostType ghost_type1;
ghost_type1 = _not_ghost;
ByElementTypeReal quadrature_points_volumes("quadrature_points_volumes", id, memory_id);
this->initInternalVector(quadrature_points_volumes, 1, true);
resizeInternalVector(quadrature_points_volumes);
const FEM & fem = this->model->getFEM();
weight_func->updateInternals(quadrature_points_volumes);
for(UInt gt = _not_ghost; gt <= _ghost; ++gt) {
GhostType ghost_type2 = (GhostType) gt;
UInt existing_pairs_num = gt - _not_ghost;
pair_type::iterator first_pair_types = existing_pairs[existing_pairs_num].begin();
pair_type::iterator last_pair_types = existing_pairs[existing_pairs_num].end();
// Compute the weights
for (; first_pair_types != last_pair_types; ++first_pair_types) {
ElementType type1 = first_pair_types->first;
ElementType type2 = first_pair_types->second;
const Vector<UInt> & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2);
std::string ghost_id = "";
if (ghost_type1 == _ghost) ghost_id = ":ghost";
ByElementTypeReal & weights_type_1 = pair_weight(type1, ghost_type1);
std::stringstream sstr; sstr << id << ":pair_weight:" << type1 << ghost_id;
weights_type_1.setID(sstr.str());
Vector<Real> * tmp_weight = NULL;
if(!weights_type_1.exists(type2, ghost_type2)) {
tmp_weight = &(weights_type_1.alloc(0, 2, type2, ghost_type2));
} else {
tmp_weight = &(weights_type_1(type2, ghost_type2));
}
Vector<Real> & weights = *tmp_weight;
weights.resize(pairs.getSize());
weights.clear();
const Vector<Real> & jacobians_1 = fem.getIntegratorInterface().getJacobians(type1, ghost_type1);
const Vector<Real> & jacobians_2 = fem.getIntegratorInterface().getJacobians(type2, ghost_type2);
UInt nb_quad1 = fem.getNbQuadraturePoints(type1);
UInt nb_quad2 = fem.getNbQuadraturePoints(type2);
Vector<Real> & quads_volumes1 = quadrature_points_volumes(type1, ghost_type1);
Vector<Real> & quads_volumes2 = quadrature_points_volumes(type2, ghost_type2);
Vector<Real>::const_iterator<types::RVector> iquads1;
Vector<Real>::const_iterator<types::RVector> iquads2;
iquads1 = quadrature_points_coordinates(type1, ghost_type1).begin(spatial_dimension);
iquads2 = quadrature_points_coordinates(type2, ghost_type2).begin(spatial_dimension);
Vector<UInt>::const_iterator< types::Vector<UInt> > first_pair = pairs.begin(2);
Vector<UInt>::const_iterator< types::Vector<UInt> > last_pair = pairs.end(2);
Vector<Real>::iterator<types::RVector> weight = weights.begin(2);
this->weight_func->selectType(type1, ghost_type1, type2, ghost_type2);
// Weight function
for(;first_pair != last_pair; ++first_pair, ++weight) {
UInt _q1 = (*first_pair)(0);
UInt _q2 = (*first_pair)(1);
const types::RVector & pos1 = iquads1[_q1];
const types::RVector & pos2 = iquads2[_q2];
QuadraturePoint q1(_q1 / nb_quad1, _q1 % nb_quad1, _q1, pos1, type1, ghost_type1);
QuadraturePoint q2(_q2 / nb_quad2, _q2 % nb_quad2, _q2, pos2, type2, ghost_type2);
Real r = pos1.distance(pos2);
Real w2J2 = jacobians_2(_q2);
(*weight)(0) = w2J2 * this->weight_func->operator()(r, q1, q2);
if(ghost_type2 != _ghost && _q1 != _q2) {
Real w1J1 = jacobians_1(_q1);
(*weight)(1) = w1J1 * this->weight_func->operator()(r, q2, q1);
} else
(*weight)(1) = 0;
quads_volumes1(_q1) += (*weight)(0);
if(ghost_type2 != _ghost) quads_volumes2(_q2) += (*weight)(1);
}
}
}
//normalize the weights
for(UInt gt = _not_ghost; gt <= _ghost; ++gt) {
GhostType ghost_type2 = (GhostType) gt;
UInt existing_pairs_num = gt - _not_ghost;
pair_type::iterator first_pair_types = existing_pairs[existing_pairs_num].begin();
pair_type::iterator last_pair_types = existing_pairs[existing_pairs_num].end();
first_pair_types = existing_pairs[existing_pairs_num].begin();
for (; first_pair_types != last_pair_types; ++first_pair_types) {
ElementType type1 = first_pair_types->first;
ElementType type2 = first_pair_types->second;
const Vector<UInt> & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2);
Vector<Real> & weights = pair_weight(type1, ghost_type1)(type2, ghost_type2);
Vector<Real> & quads_volumes1 = quadrature_points_volumes(type1, ghost_type1);
Vector<Real> & quads_volumes2 = quadrature_points_volumes(type2, ghost_type2);
Vector<UInt>::const_iterator< types::Vector<UInt> > first_pair = pairs.begin(2);
Vector<UInt>::const_iterator< types::Vector<UInt> > last_pair = pairs.end(2);
Vector<Real>::iterator<types::RVector> weight = weights.begin(2);
for(;first_pair != last_pair; ++first_pair, ++weight) {
UInt q1 = (*first_pair)(0);
UInt q2 = (*first_pair)(1);
(*weight)(0) *= 1. / quads_volumes1(q1);
if(ghost_type2 != _ghost) (*weight)(1) *= 1. / quads_volumes2(q2);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
template<typename T>
void MaterialNonLocal<spatial_dimension, WeightFunction>::weightedAvergageOnNeighbours(const ByElementTypeVector<T> & to_accumulate,
ByElementTypeVector<T> & accumulated,
UInt nb_degree_of_freedom,
GhostType ghost_type2) const {
AKANTU_DEBUG_IN();
UInt existing_pairs_num = 0;
if (ghost_type2 == _ghost) existing_pairs_num = 1;
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();
GhostType ghost_type1 = _not_ghost; // does not make sens the ghost vs ghost so this should always by _not_ghost
for (; first_pair_types != last_pair_types; ++first_pair_types) {
const Vector<UInt> & pairs =
pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2);
const Vector<Real> & weights =
pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2);
const Vector<T> & to_acc = to_accumulate(first_pair_types->second, ghost_type2);
Vector<T> & acc = accumulated(first_pair_types->first, ghost_type1);
if(ghost_type2 == _not_ghost) acc.clear();
Vector<UInt>::const_iterator< types::Vector<UInt> > first_pair = pairs.begin(2);
Vector<UInt>::const_iterator< types::Vector<UInt> > last_pair = pairs.end(2);
Vector<Real>::const_iterator< types::Vector<Real> > pair_w = weights.begin(2);
for(;first_pair != last_pair; ++first_pair, ++pair_w) {
UInt q1 = (*first_pair)(0);
UInt q2 = (*first_pair)(1);
for(UInt d = 0; d < nb_degree_of_freedom; ++d){
acc(q1, d) += (*pair_w)(0) * to_acc(q2, d);
if(ghost_type2 != _ghost) acc(q2, d) += (*pair_w)(1) * to_acc(q1, d);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::updateResidual(GhostType ghost_type) {
AKANTU_DEBUG_IN();
// Update the weights for the non local variable averaging
if(ghost_type == _not_ghost &&
this->update_weights &&
(this->compute_stress_calls % this->update_weights == 0)) {
ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id);
Mesh & mesh = this->model->getFEM().getMesh();
mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0);
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<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::computeAllNonLocalStresses(GhostType ghost_type) {
// Update the weights for the non local variable averaging
if(ghost_type == _not_ghost) {
if(this->update_weights && (this->compute_stress_calls % this->update_weights == 0)) {
this->model->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_weight);
ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id);
Mesh & mesh = this->model->getFEM().getMesh();
mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0);
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<ID, NonLocalVariable>::iterator it = non_local_variables.begin();
typename std::map<ID, NonLocalVariable>::iterator end = non_local_variables.end();
for(;it != end; ++it) {
NonLocalVariable & non_local_variable = it->second;
resizeInternalVector(*non_local_variable.non_local_variable);
this->weightedAvergageOnNeighbours(*non_local_variable.local_variable, *non_local_variable.non_local_variable,
non_local_variable.non_local_variable_nb_component, _not_ghost);
}
++this->compute_stress_calls;
} else {
typename std::map<ID, NonLocalVariable>::iterator it = non_local_variables.begin();
typename std::map<ID, NonLocalVariable>::iterator end = non_local_variables.end();
for(;it != end; ++it) {
NonLocalVariable & non_local_variable = it->second;
this->weightedAvergageOnNeighbours(*non_local_variable.local_variable, *non_local_variable.non_local_variable,
non_local_variable.non_local_variable_nb_component, _ghost);
}
computeNonLocalStresses(_not_ghost);
}
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
bool MaterialNonLocal<spatial_dimension, WeightFunction>::parseParam(const std::string & key,
const std::string & value,
__attribute__((unused)) const ID & id) {
std::stringstream sstr(value);
if(key == "radius") { sstr >> radius; }
else if(key == "UpdateWeights") { sstr >> update_weights; }
else if(!weight_func->parseParam(key, value)) return Material::parseParam(key, value, id);
return true;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::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());
GhostType ghost_type1;
ghost_type1 = _not_ghost;
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 Vector<UInt> & pairs =
pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2);
const Vector<Real> & weights =
pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2);
pout << "Types : " << first_pair_types->first << " (" << ghost_type1 << ") - " << first_pair_types->second << " (" << ghost_type2 << ")" << std::endl;
Vector<UInt>::const_iterator< types::Vector<UInt> > first_pair = pairs.begin(2);
Vector<UInt>::const_iterator< types::Vector<UInt> > last_pair = pairs.end(2);
Vector<Real>::const_iterator<types::RVector> pair_w = weights.begin(2);
for(;first_pair != last_pair; ++first_pair, ++pair_w) {
UInt q1 = (*first_pair)(0);
UInt q2 = (*first_pair)(1);
pout << q1 << " " << q2 << " "<< (*pair_w)(0) << " " << (*pair_w)(1) << std::endl;
}
}
}
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
void MaterialNonLocal<spatial_dimension, WeightFunction>::neighbourhoodStatistics(const std::string & filename) const {
std::ofstream pout;
pout.open(filename.c_str());
const Mesh & mesh = this->model->getFEM().getMesh();
GhostType ghost_type1;
ghost_type1 = _not_ghost;
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
Int prank = comm.whoAmI();
for(UInt gt = _not_ghost; gt <= _ghost; ++gt) {
GhostType ghost_type2 = (GhostType) gt;
UInt existing_pairs_num = gt - _not_ghost;
ByElementTypeUInt nb_neighbors("nb_neighbours", id, memory_id);
mesh.initByElementTypeVector(nb_neighbors, 1, spatial_dimension);
resizeInternalVector(nb_neighbors);
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 Vector<UInt> & 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;
}
Vector<UInt>::const_iterator< types::Vector<UInt> > first_pair = pairs.begin(2);
Vector<UInt>::const_iterator< types::Vector<UInt> > last_pair = pairs.end(2);
Vector<UInt> & nb_neigh_1 = nb_neighbors(first_pair_types->first, ghost_type1);
Vector<UInt> & 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<UInt>::max();
for(; it != last_type; ++it) {
Vector<UInt> & nb_neighor = nb_neighbors(*it, ghost_type1);
Vector<UInt>::iterator<UInt> nb_neigh = nb_neighor.begin();
Vector<UInt>::iterator<UInt> 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<UInt spatial_dimension, template <UInt> class WeightFunction>
inline UInt MaterialNonLocal<spatial_dimension, WeightFunction>::getNbDataForElements(const Vector<Element> & elements,
SynchronizationTag tag) const {
UInt nb_quadrature_points = this->getModel().getNbQuadraturePoints(elements);
UInt size = 0;
if(tag == _gst_mnl_for_average) {
typename std::map<ID, NonLocalVariable>::const_iterator it = non_local_variables.begin();
typename std::map<ID, NonLocalVariable>::const_iterator end = non_local_variables.end();
for(;it != end; ++it) {
const NonLocalVariable & non_local_variable = it->second;
size += non_local_variable.non_local_variable_nb_component * sizeof(Real) * nb_quadrature_points;
}
}
size += weight_func->getNbDataForElements(elements, tag);
return size;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
inline void MaterialNonLocal<spatial_dimension, WeightFunction>::packElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag) const {
if(tag == _gst_mnl_for_average) {
typename std::map<ID, NonLocalVariable>::const_iterator it = non_local_variables.begin();
typename std::map<ID, NonLocalVariable>::const_iterator end = non_local_variables.end();
for(;it != end; ++it) {
const NonLocalVariable & non_local_variable = it->second;
this->packElementDataHelper(*non_local_variable.local_variable,
buffer, elements);
}
}
weight_func->packElementData(buffer, elements, tag);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
inline void MaterialNonLocal<spatial_dimension, WeightFunction>::unpackElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag) {
if(tag == _gst_mnl_for_average) {
typename std::map<ID, NonLocalVariable>::iterator it = non_local_variables.begin();
typename std::map<ID, NonLocalVariable>::iterator end = non_local_variables.end();
for(;it != end; ++it) {
NonLocalVariable & non_local_variable = it->second;
this->unpackElementDataHelper(*non_local_variable.local_variable,
buffer, elements);
}
}
weight_func->unpackElementData(buffer, elements, tag);
}
/* -------------------------------------------------------------------------- */
// template<UInt spatial_dimension, template <UInt> class WeightFunction>
// inline void MaterialNonLocal<spatial_dimension, WeightFunction>::onElementsAdded(const Vector<Element> & element_list) {
// AKANTU_DEBUG_IN();
// Material::onElementsAdded(element_list, event);
// AKANTU_DEBUG_OUT();
// }
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension, template <UInt> class WeightFunction>
inline void MaterialNonLocal<spatial_dimension, WeightFunction>::onElementsRemoved(const Vector<Element> & element_list,
const ByElementTypeUInt & new_numbering,
__attribute__((unused)) const RemovedElementsEvent & event) {
AKANTU_DEBUG_IN();
Material::onElementsRemoved(element_list, new_numbering, event);
pair_type::const_iterator first_pair_types = existing_pairs[1].begin();
pair_type::const_iterator last_pair_types = existing_pairs[1].end();
// Renumber element to keep
for (; first_pair_types != last_pair_types; ++first_pair_types) {
ElementType type2 = first_pair_types->second;
GhostType ghost_type2 = _ghost;
UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2);
Vector<UInt> & pairs =
pair_list(first_pair_types->first, _not_ghost)(first_pair_types->second, ghost_type2);
Vector<UInt>::iterator< types::Vector<UInt> > first_pair = pairs.begin(2);
Vector<UInt>::iterator< types::Vector<UInt> > last_pair = pairs.end(2);
for(;first_pair != last_pair; ++first_pair) {
UInt _q2 = (*first_pair)(1);
const Vector<UInt> & renumbering = new_numbering(type2, ghost_type2);
UInt el = _q2 / nb_quad2;
UInt new_el = renumbering(el);
AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A local quad as been removed instead f just renumbered");
(*first_pair)(1) = new_el * nb_quad2 + _q2 % nb_quad2;
}
}
AKANTU_DEBUG_OUT();
}
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index 1c9f71b17..79b9b615e 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1390 +1,1394 @@
/**
* @file solid_mechanics_model.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Tue Jul 27 18:15:37 2010
*
* @brief Implementation of the SolidMechanicsModel class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_math.hh"
#include "aka_common.hh"
#include "solid_mechanics_model.hh"
#include "integration_scheme_2nd_order.hh"
#include "static_communicator.hh"
#include "dof_synchronizer.hh"
#include <cmath>
#ifdef AKANTU_USE_MUMPS
#include "solver_mumps.hh"
#endif
-/* -------------------------------------------------------------------------- */
-
+/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
+#ifdef AKANTU_USE_IOHELPER
+# include "dumper_iohelper_tmpl_homogenizing_field.hh"
+# include "dumper_iohelper_tmpl_material_internal_field.hh"
+#endif
/* -------------------------------------------------------------------------- */
/**
* A solid mechanics model need a mesh and a dimension to be created. the model
* by it self can not do a lot, the good init functions should be called in
* order to configure the model depending on what we want to do.
*
* @param mesh mesh representing the model we want to simulate
* @param dim spatial dimension of the problem, if dim = 0 (default value) the
* dimension of the problem is assumed to be the on of the mesh
* @param id an id to identify the model
*/
SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh,
UInt dim,
const ID & id,
const MemoryID & memory_id) :
Model(mesh, id, memory_id), Dumpable<DumperParaview>(id),
time_step(NAN), f_m2a(1.0),
mass_matrix(NULL),
velocity_damping_matrix(NULL),
stiffness_matrix(NULL),
jacobian_matrix(NULL),
//element_material("element_material", id),
element_index_by_material("element index by material", id),
integrator(NULL),
increment_flag(false), solver(NULL),
spatial_dimension(dim) {
AKANTU_DEBUG_IN();
createSynchronizerRegistry(this);
if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension();
registerFEMObject<MyFEMType>("SolidMechanicsFEM", mesh, spatial_dimension);
this->displacement = NULL;
this->mass = NULL;
this->velocity = NULL;
this->acceleration = NULL;
this->force = NULL;
this->residual = NULL;
this->boundary = NULL;
this->increment = NULL;
this->increment_acceleration = NULL;
this->dof_synchronizer = NULL;
materials.clear();
mesh.registerEventHandler(*this);
addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::~SolidMechanicsModel() {
AKANTU_DEBUG_IN();
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
delete *mat_it;
}
materials.clear();
delete integrator;
if(solver) delete solver;
if(mass_matrix) delete mass_matrix;
if(velocity_damping_matrix) delete velocity_damping_matrix;
if(stiffness_matrix && stiffness_matrix != jacobian_matrix)
delete stiffness_matrix;
if(jacobian_matrix) delete jacobian_matrix;
if(dof_synchronizer) delete dof_synchronizer;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Initialisation */
/* -------------------------------------------------------------------------- */
/**
* This function groups many of the initialization in on function. For most of
* basics case the function should be enough. The functions initialize the
* model, the internal vectors, set them to 0, and depending on the parameters
* it also initialize the explicit or implicit solver.
*
* @param material_file the file containing the materials to use
* @param method the analysis method wanted. See the akantu::AnalysisMethod for
* the different possibilites
*/
void SolidMechanicsModel::initFull(std::string material_file,
AnalysisMethod analysis_method) {
method = analysis_method;
// initialize the model
initModel();
// initialize the vectors
initVectors();
// set the initial condition to 0
force->clear();
velocity->clear();
acceleration->clear();
displacement->clear();
// initialize pcb
if(pbc_pair.size()!=0)
initPBC();
// initialize the time integration schemes
switch(method) {
case _explicit_dynamic:
initExplicit();
break;
case _implicit_dynamic:
initImplicit(true);
break;
case _static:
initImplicit(false);
break;
}
// initialize the materials
if(material_file != "") {
readMaterials(material_file);
initMaterials();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initParallel(MeshPartition * partition,
DataAccessor * data_accessor) {
AKANTU_DEBUG_IN();
if (data_accessor == NULL) data_accessor = this;
synch_parallel = &createParallelSynch(partition,data_accessor);
synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass);
// synch_registry->registerSynchronizer(synch_parallel, _gst_smm_for_strain);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initFEMBoundary(bool create_surface) {
if(create_surface)
MeshUtils::buildFacets(mesh);
FEM & fem_boundary = getFEMBoundary();
fem_boundary.initShapeFunctions();
fem_boundary.computeNormalsOnControlPoints();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initExplicit() {
AKANTU_DEBUG_IN();
method = _explicit_dynamic;
if (integrator) delete integrator;
integrator = new CentralDifference();
UInt nb_nodes = acceleration->getSize();
UInt nb_degree_of_freedom = acceleration->getNbComponent();
std::stringstream sstr; sstr << id << ":increment_acceleration";
increment_acceleration = &(alloc<Real>(sstr.str(), nb_nodes, nb_degree_of_freedom, Real()));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Allocate all the needed vectors. By default their are not necessarily set to
* 0
*
*/
void SolidMechanicsModel::initVectors() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_disp; sstr_disp << id << ":displacement";
std::stringstream sstr_mass; sstr_mass << id << ":mass";
std::stringstream sstr_velo; sstr_velo << id << ":velocity";
std::stringstream sstr_acce; sstr_acce << id << ":acceleration";
std::stringstream sstr_forc; sstr_forc << id << ":force";
std::stringstream sstr_resi; sstr_resi << id << ":residual";
std::stringstream sstr_boun; sstr_boun << id << ":boundary";
displacement = &(alloc<Real>(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
mass = &(alloc<Real>(sstr_mass.str(), nb_nodes, spatial_dimension, 0));
velocity = &(alloc<Real>(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
acceleration = &(alloc<Real>(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
force = &(alloc<Real>(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
residual = &(alloc<Real>(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
boundary = &(alloc<bool>(sstr_boun.str(), nb_nodes, spatial_dimension, false));
std::stringstream sstr_curp; sstr_curp << id << ":current_position";
current_position = &(alloc<Real>(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE));
for(UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType) g;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined);
Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined);
for(; it != end; ++it) {
UInt nb_element = mesh.getNbElement(*it, gt);
element_index_by_material.alloc(nb_element, 2, *it, gt);
}
}
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*
*/
void SolidMechanicsModel::initModel() {
/// \todo add the current position as a parameter to initShapeFunctions for
/// large deformation
getFEM().initShapeFunctions(_not_ghost);
getFEM().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initPBC() {
Model::initPBC();
registerPBCSynchronizer();
// as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves
std::map<UInt, UInt>::iterator it = pbc_pair.begin();
std::map<UInt, UInt>::iterator end = pbc_pair.end();
UInt dim = mesh.getSpatialDimension();
while(it != end) {
for (UInt i=0; i<dim; ++i)
(*boundary)((*it).first,i) = true;
++it;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::registerPBCSynchronizer(){
PBCSynchronizer * synch = new PBCSynchronizer(pbc_pair);
synch_registry->registerSynchronizer(*synch, _gst_smm_uv);
synch_registry->registerSynchronizer(*synch, _gst_smm_mass);
synch_registry->registerSynchronizer(*synch, _gst_smm_res);
changeLocalEquationNumberforPBC(pbc_pair, mesh.getSpatialDimension());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateCurrentPosition() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
current_position->resize(nb_nodes);
Real * current_position_val = current_position->values;
Real * position_val = mesh.getNodes().values;
Real * displacement_val = displacement->values;
/// compute current_position = initial_position + displacement
memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real));
for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) {
*current_position_val++ += *displacement_val++;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initializeUpdateResidualData() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
residual->resize(nb_nodes);
/// copy the forces in residual for boundary conditions
memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real));
// start synchronization
synch_registry->asynchronousSynchronize(_gst_smm_uv);
synch_registry->waitEndSynchronize(_gst_smm_uv);
updateCurrentPosition();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Explicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* This function compute the second member of the motion equation. That is to
* say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given
* by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$
* F_{int} = \int_{\Omega} N \sigma d\Omega@f$
*
*/
void SolidMechanicsModel::updateResidual(bool need_initialize) {
AKANTU_DEBUG_IN();
// f = f_ext - f_int
// f = f_ext
if (need_initialize) initializeUpdateResidualData();
// communicate the displacement
// synch_registry->asynchronousSynchronize(_gst_smm_for_strain);
std::vector<Material *>::iterator mat_it;
// call update residual on each local elements
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStresses(_not_ghost);
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* ------------------------------------------------------------------------ */
/* Computation of the non local part */
synch_registry->asynchronousSynchronize(_gst_mnl_for_average);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_not_ghost);
}
synch_registry->waitEndSynchronize(_gst_mnl_for_average);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_ghost);
}
#endif
/* ------------------------------------------------------------------------ */
/* assembling the forces internal */
// communicate the strain
synch_registry->asynchronousSynchronize(_gst_smm_stress);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleResidual(_not_ghost);
}
// finalize communications
synch_registry->waitEndSynchronize(_gst_smm_stress);
// synch_registry->waitEndSynchronize(_gst_smm_for_strain);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleResidual(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeStresses() {
if (method == _explicit_dynamic) {
// start synchronization
synch_registry->asynchronousSynchronize(_gst_smm_uv);
synch_registry->waitEndSynchronize(_gst_smm_uv);
// compute stresses on all local elements for each materials
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStresses(_not_ghost);
}
/* ------------------------------------------------------------------------ */
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* Computation of the non local part */
synch_registry->asynchronousSynchronize(_gst_mnl_for_average);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_not_ghost);
}
synch_registry->waitEndSynchronize(_gst_mnl_for_average);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_ghost);
}
#endif
} else {
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStressesFromTangentModuli(_not_ghost);
}
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateResidualInternal() {
AKANTU_DEBUG_IN();
// f = f_ext - f_int - Ma - Cv = r - Ma - Cv;
if(method != _static) {
// f -= Ma
if(mass_matrix) {
// if full mass_matrix
Vector<Real> * Ma = new Vector<Real>(*acceleration, true, "Ma");
*Ma *= *mass_matrix;
/// \todo check unit conversion for implicit dynamics
// *Ma /= f_m2a
*residual -= *Ma;
delete Ma;
} else {
// else lumped mass
UInt nb_nodes = acceleration->getSize();
UInt nb_degree_of_freedom = acceleration->getNbComponent();
Real * mass_val = mass->values;
Real * accel_val = acceleration->values;
Real * res_val = residual->values;
bool * boundary_val = boundary->values;
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if(!(*boundary_val)) {
*res_val -= *accel_val * *mass_val /f_m2a;
}
boundary_val++;
res_val++;
mass_val++;
accel_val++;
}
}
// f -= Cv
if(velocity_damping_matrix) {
Vector<Real> * Cv = new Vector<Real>(*velocity);
*Cv *= *velocity_damping_matrix;
*residual -= *Cv;
delete Cv;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateAcceleration() {
AKANTU_DEBUG_IN();
updateResidualInternal();
if(method == _explicit_dynamic && !mass_matrix) {
/* residual = residual_{n+1} - M * acceleration_n therefore
solution = increment acceleration not acceleration */
solveLumped(*increment_acceleration,
*mass,
*residual,
*boundary,
f_m2a);
} else {
solveDynamic<NewmarkBeta::_acceleration_corrector>(*increment_acceleration);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveLumped(Vector<Real> & x,
const Vector<Real> & A,
const Vector<Real> & b,
const Vector<bool> & boundary,
Real alpha) {
Real * A_val = A.storage();
Real * b_val = b.storage();
Real * x_val = x.storage();
bool * boundary_val = boundary.storage();
UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent();
for (UInt n = 0; n < nb_degrees_of_freedom; ++n) {
if(!(*boundary_val)) {
*x_val = alpha * (*b_val / *A_val);
}
x_val++;
A_val++;
b_val++;
boundary_val++;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitPred() {
AKANTU_DEBUG_IN();
if(increment_flag) {
memcpy(increment->values,
displacement->values,
displacement->getSize()*displacement->getNbComponent()*sizeof(Real));
}
AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: "
<< "have called initExplicit ? "
<< "or initImplicit ?");
integrator->integrationSchemePred(time_step,
*displacement,
*velocity,
*acceleration,
*boundary);
if(increment_flag) {
Real * inc_val = increment->values;
Real * dis_val = displacement->values;
UInt nb_nodes = displacement->getSize();
for (UInt n = 0; n < nb_nodes; ++n) {
*inc_val = *dis_val - *inc_val;
inc_val++;
dis_val++;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitCorr() {
AKANTU_DEBUG_IN();
integrator->integrationSchemeCorrAccel(time_step,
*displacement,
*velocity,
*acceleration,
*boundary,
*increment_acceleration);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Implicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* Initialize the solver and create the sparse matrices needed.
*
*/
void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) {
#if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#else
UInt nb_global_node = mesh.getNbGlobalNodes();
std::stringstream sstr; sstr << id << ":jacobian_matrix";
jacobian_matrix = new SparseMatrix(nb_global_node * spatial_dimension, _symmetric,
spatial_dimension, sstr.str(), memory_id);
jacobian_matrix->buildProfile(mesh, *dof_synchronizer);
std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix";
stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id);
#ifdef AKANTU_USE_MUMPS
std::stringstream sstr_solv; sstr_solv << id << ":solver";
solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());
dof_synchronizer->initScatterGatherCommunicationScheme();
#else
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#endif //AKANTU_USE_MUMPS
solver->initialize(options);
#endif //AKANTU_HAS_SOLVER
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the implicit solver, either for dynamic or static cases,
*
* @param dynamic
*/
void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) {
AKANTU_DEBUG_IN();
method = dynamic ? _implicit_dynamic : _static;
initSolver(solver_options);
if(method == _implicit_dynamic) {
if(integrator) delete integrator;
integrator = new TrapezoidalRule2();
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initialAcceleration() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Solving Ma = f");
Solver * acc_solver = NULL;
std::stringstream sstr; sstr << id << ":tmp_mass_matrix";
SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id);
#ifdef AKANTU_USE_MUMPS
std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix";
acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str());
dof_synchronizer->initScatterGatherCommunicationScheme();
#else
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#endif //AKANTU_USE_MUMPS
acc_solver->initialize();
tmp_mass->applyBoundary(*boundary);
acc_solver->setRHS(*residual);
acc_solver->solve(*acceleration);
delete acc_solver;
delete tmp_mass;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
stiffness_matrix->clear();
// call compute stiffness matrix on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->assembleStiffnessMatrix(_not_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveDynamic() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL,
"You should first initialize the implicit solver and assemble the stiffness matrix");
AKANTU_DEBUG_ASSERT(mass_matrix != NULL,
"You should first initialize the implicit solver and assemble the mass matrix");
if(!increment) setIncrementFlagOn();
updateResidualInternal();
solveDynamic<NewmarkBeta::_displacement_corrector>(*increment);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveStatic() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Solving Ku = f");
AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL,
"You should first initialize the implicit solver and assemble the stiffness matrix");
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent();
// if(method != _static)
jacobian_matrix->copyContent(*stiffness_matrix);
jacobian_matrix->applyBoundary(*boundary);
solver->setRHS(*residual);
if(!increment) setIncrementFlagOn();
solver->solve(*increment);
Real * increment_val = increment->values;
Real * displacement_val = displacement->values;
bool * boundary_val = boundary->values;
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if(!(*boundary_val)) {
*displacement_val += *increment_val;
}
displacement_val++;
boundary_val++;
increment_val++;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance) {
Real error;
bool tmp = testConvergenceIncrement(tolerance, error);
AKANTU_DEBUG_INFO("Norm of increment : " << error);
return tmp;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error) {
AKANTU_DEBUG_IN();
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent();
error = 0;
Real norm[2] = {0., 0.};
Real * increment_val = increment->storage();
bool * boundary_val = boundary->storage();
Real * displacement_val = displacement->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
if(!(*boundary_val) && is_local_node) {
norm[0] += *increment_val * *increment_val;
norm[1] += *displacement_val * *displacement_val;
}
boundary_val++;
increment_val++;
displacement_val++;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum);
norm[0] = sqrt(norm[0]);
norm[1] = sqrt(norm[1]);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
error = norm[0] / norm[1];
return (error < tolerance);
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceResidual(Real tolerance) {
Real error;
bool tmp = testConvergenceResidual(tolerance, error);
AKANTU_DEBUG_INFO("Norm of residual : " << error);
return tmp;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & norm) {
AKANTU_DEBUG_IN();
UInt nb_nodes = residual->getSize();
norm = 0;
Real * residual_val = residual->values;
bool * boundary_val = boundary->values;
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
if(is_local_node) {
for (UInt d = 0; d < spatial_dimension; ++d) {
if(!(*boundary_val)) {
norm += *residual_val * *residual_val;
}
boundary_val++;
residual_val++;
}
} else {
boundary_val += spatial_dimension;
residual_val += spatial_dimension;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum);
norm = sqrt(norm);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
return (norm < tolerance);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::implicitPred() {
AKANTU_DEBUG_IN();
integrator->integrationSchemePred(time_step,
*displacement,
*velocity,
*acceleration,
*boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::implicitCorr() {
AKANTU_DEBUG_IN();
integrator->integrationSchemeCorrDispl(time_step,
*displacement,
*velocity,
*acceleration,
*boundary,
*increment);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Information */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeBoundaries() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized."
<< " Did you call initParallel?");
synch_registry->synchronize(_gst_smm_boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeResidual() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized."
<< " Did you call initPBC?");
synch_registry->synchronize(_gst_smm_res);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setIncrementFlagOn() {
AKANTU_DEBUG_IN();
if(!increment) {
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_inc; sstr_inc << id << ":increment";
increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, spatial_dimension, 0));
}
increment_flag = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep() {
AKANTU_DEBUG_IN();
Real min_dt = getStableTimeStep(_not_ghost);
/// reduction min over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
Material ** mat_val = &(materials.at(0));
Real min_dt = std::numeric_limits<Real>::max();
updateCurrentPosition();
Element elem;
elem.ghost_type = ghost_type;
elem.kind = _ek_regular;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type);
for(; it != end; ++it) {
elem.type = *it;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it);
UInt nb_element = mesh.getNbElement(*it);
Vector<UInt>::iterator< types::Vector<UInt> > eibm =
element_index_by_material(*it, ghost_type).begin(2);
Vector<Real> X(0, nb_nodes_per_element*spatial_dimension);
FEM::extractNodalToElementField(mesh, *current_position,
X, *it, _not_ghost);
Vector<Real>::iterator< types::Matrix<Real> > X_el =
X.begin(spatial_dimension, nb_nodes_per_element);
for (UInt el = 0; el < nb_element; ++el, ++X_el, ++eibm) {
elem.element = el;
Real el_size = getFEM().getElementInradius(*X_el, *it);
Real el_dt = mat_val[(*eibm)(1)]->getStableTimeStep(el_size, elem);
min_dt = std::min(min_dt, el_dt);
}
}
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getPotentialEnergy() {
AKANTU_DEBUG_IN();
Real energy = 0.;
/// call update residual on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
energy += (*mat_it)->getPotentialEnergy();
}
/// reduction sum over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy() {
AKANTU_DEBUG_IN();
Real ekin = 0.;
UInt nb_nodes = mesh.getNbNodes();
Real * vel_val = velocity->values;
Real * mass_val = mass->values;
for (UInt n = 0; n < nb_nodes; ++n) {
Real mv2 = 0;
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < spatial_dimension; ++i) {
if (count_node)
mv2 += *vel_val * *vel_val * *mass_val;
vel_val++;
mass_val++;
}
ekin += mv2;
}
StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum);
AKANTU_DEBUG_OUT();
return ekin * .5;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getExternalWork() {
AKANTU_DEBUG_IN();
Real * velo = velocity->storage();
Real * forc = force->storage();
Real * resi = residual->storage();
bool * boun = boundary->storage();
Real work = 0.;
UInt nb_nodes = mesh.getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < spatial_dimension; ++i) {
if (count_node) {
if(*boun)
work -= *resi * *velo * time_step;
else
work += *forc * *velo * time_step;
}
++velo;
++forc;
++resi;
++boun;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum);
AKANTU_DEBUG_OUT();
return work;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & id) {
AKANTU_DEBUG_IN();
if (id == "kinetic") {
return getKineticEnergy();
} else if (id == "external work"){
return getExternalWork();
}
Real energy = 0.;
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
energy += (*mat_it)->getEnergy(id);
}
/// reduction sum over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
ElementType & type,
UInt index){
AKANTU_DEBUG_IN();
std::vector<Material *>::iterator mat_it;
types::Vector<UInt> mat = element_index_by_material(type, _not_ghost).begin(2)[index];
Real energy = materials[mat(1)]->getEnergy(energy_id, type, mat(0));
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesAdded(const Vector<UInt> & nodes_list,
__attribute__((unused)) const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
if(displacement) displacement->resize(nb_nodes);
if(mass ) mass ->resize(nb_nodes);
if(velocity ) velocity ->resize(nb_nodes);
if(acceleration) acceleration->resize(nb_nodes);
if(force ) force ->resize(nb_nodes);
if(residual ) residual ->resize(nb_nodes);
if(boundary ) boundary ->resize(nb_nodes);
if(increment_acceleration) increment_acceleration->resize(nb_nodes);
if(increment) increment->resize(nb_nodes);
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onNodesAdded(nodes_list, event);
}
if (method != _explicit_dynamic) {
delete stiffness_matrix;
delete jacobian_matrix;
delete solver;
SolverOptions solver_options;
initImplicit((method == _implicit_dynamic), solver_options);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsAdded(const Vector<Element> & element_list,
const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
getFEM().initShapeFunctions(_not_ghost);
getFEM().initShapeFunctions(_ghost);
Vector<Element>::const_iterator<Element> it = element_list.begin();
Vector<Element>::const_iterator<Element> end = element_list.end();
/// \todo have rules to choose the correct material
UInt mat_id = 0;
UInt * mat_id_vect = NULL;
try {
const NewMaterialElementsEvent & event_mat = dynamic_cast<const NewMaterialElementsEvent &>(event);
mat_id_vect = event_mat.getMaterialList().storage();
} catch(...) { }
for (UInt el = 0; it != end; ++it, ++el) {
const Element & elem = *it;
// element_material(elem.type, elem.ghost_type).push_back(UInt(0));
if(mat_id_vect) mat_id = mat_id_vect[el];
Material & mat = *materials[mat_id];
UInt mat_index = mat.addElement(elem.type, elem.element, elem.ghost_type);
UInt id[2];
id[0] = mat_index; id[1] = 0;
element_index_by_material(elem.type, elem.ghost_type).push_back(id);
}
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsAdded(element_list, event);
}
if(method != _explicit_dynamic) AKANTU_DEBUG_TO_IMPLEMENT();
assembleMassLumped();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Vector<Element> & element_list,
const ByElementTypeUInt & new_numbering,
const RemovedElementsEvent & event) {
// MeshUtils::purifyMesh(mesh);
getFEM().initShapeFunctions(_not_ghost);
getFEM().initShapeFunctions(_ghost);
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsRemoved(element_list, new_numbering, event);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Vector<UInt> & element_list,
const Vector<UInt> & new_numbering,
__attribute__((unused)) const RemovedNodesEvent & event) {
if(displacement) mesh.removeNodesFromVector(*displacement, new_numbering);
if(mass ) mesh.removeNodesFromVector(*mass , new_numbering);
if(velocity ) mesh.removeNodesFromVector(*velocity , new_numbering);
if(acceleration) mesh.removeNodesFromVector(*acceleration, new_numbering);
if(force ) mesh.removeNodesFromVector(*force , new_numbering);
if(residual ) mesh.removeNodesFromVector(*residual , new_numbering);
if(boundary ) mesh.removeNodesFromVector(*boundary , new_numbering);
if(increment_acceleration) mesh.removeNodesFromVector(*increment_acceleration, new_numbering);
if(increment) mesh.removeNodesFromVector(*increment , new_numbering);
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpField(const std::string & field_id) {
#ifdef AKANTU_USE_IOHELPER
#define ADD_FIELD(field, type) \
addDumpFieldToDumper(BOOST_PP_STRINGIZE(field), \
new DumperIOHelper::NodalField<type>(*field))
if(field_id == "displacement") { ADD_FIELD(displacement, Real); }
else if(field_id == "mass" ) { ADD_FIELD(mass , Real); }
else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); }
else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); }
else if(field_id == "force" ) { ADD_FIELD(force , Real); }
else if(field_id == "residual" ) { ADD_FIELD(residual , Real); }
else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); }
else if(field_id == "partitions" ) {
addDumpFieldToDumper(field_id,
new DumperIOHelper::ElementPartitionField(mesh,
spatial_dimension,
_not_ghost,
_ek_regular));
}
else if(field_id == "element_index_by_material") {
addDumpFieldToDumper(field_id,
new DumperIOHelper::ElementalField<UInt>(element_index_by_material,
spatial_dimension,
_not_ghost,
_ek_regular));
} else {
addDumpFieldToDumper(field_id,
new DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField>(*this,
field_id,
spatial_dimension,
_not_ghost,
_ek_regular));
}
#undef ADD_FIELD
#endif
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpFieldVector(const std::string & field_id) {
#ifdef AKANTU_USE_IOHELPER
#define ADD_FIELD(field, type) \
DumperIOHelper::Field * f = \
new DumperIOHelper::NodalField<type>(*field); \
f->setPadding(3); \
addDumpFieldToDumper(BOOST_PP_STRINGIZE(field), f)
if(field_id == "displacement") { ADD_FIELD(displacement, Real); }
else if(field_id == "mass" ) { ADD_FIELD(mass , Real); }
else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); }
else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); }
else if(field_id == "force" ) { ADD_FIELD(force , Real); }
else if(field_id == "residual" ) { ADD_FIELD(residual , Real); }
else {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField> Field;
Field * field = new Field(*this,
field_id,
spatial_dimension,
spatial_dimension,
_not_ghost,
_ek_regular);
field->setPadding(3);
addDumpFieldToDumper(field_id, field);
}
#undef ADD_FIELD
#endif
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpFieldTensor(const std::string & field_id) {
#ifdef AKANTU_USE_IOHELPER
if(field_id == "stress") {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField,
+ DumperIOHelper::material_stress_field_iterator,
DumperIOHelper::AvgHomogenizingFunctor,
- types::Matrix,
- DumperIOHelper::material_stress_field_iterator> Field;
+ types::Matrix> Field;
Field * field = new Field(*this,
field_id,
spatial_dimension,
spatial_dimension,
_not_ghost,
_ek_regular);
field->setPadding(3, 3);
addDumpFieldToDumper(field_id, field);
} else if(field_id == "strain") {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField,
+ DumperIOHelper::material_strain_field_iterator,
DumperIOHelper::AvgHomogenizingFunctor,
- types::Matrix,
- DumperIOHelper::material_strain_field_iterator> Field;
+ types::Matrix> Field;
Field * field = new Field(*this,
field_id,
spatial_dimension,
spatial_dimension,
_not_ghost,
_ek_regular);
field->setPadding(3, 3);
addDumpFieldToDumper(field_id, field);
} else {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField,
+ DumperIOHelper::internal_material_field_iterator,
DumperIOHelper::AvgHomogenizingFunctor,
types::Matrix> Field;
Field * field = new Field(*this,
field_id,
spatial_dimension,
spatial_dimension,
_not_ghost,
_ek_regular);
field->setPadding(3, 3);
addDumpFieldToDumper(field_id, field);
}
#endif
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "Solid Mechanics Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << spatial_dimension << std::endl;
stream << space << " + fem [" << std::endl;
getFEM().printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
mass ->printself(stream, indent + 2);
velocity ->printself(stream, indent + 2);
acceleration->printself(stream, indent + 2);
force ->printself(stream, indent + 2);
residual ->printself(stream, indent + 2);
boundary ->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + connectivity type information [" << std::endl;
element_index_by_material.printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh
index 84651689d..c7a506e0d 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model.hh
@@ -1,594 +1,588 @@
/**
* @file solid_mechanics_model.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Tue Jul 27 18:15:37 2010
*
* @brief Model of Solid Mechanics
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__
#define __AKANTU_SOLID_MECHANICS_MODEL_HH__
/* -------------------------------------------------------------------------- */
#include <fstream>
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "model.hh"
#include "data_accessor.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
#include "aka_types.hh"
#include "integration_scheme_2nd_order.hh"
#include "solver.hh"
#include "dumpable.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
class Material;
class IntegrationScheme2ndOrder;
class Contact;
class SparseMatrix;
}
__BEGIN_AKANTU__
class SolidMechanicsModel : public Model, public DataAccessor, public MeshEventHandler, public Dumpable<DumperParaview> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
class NewMaterialElementsEvent : public NewElementsEvent {
public:
AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Vector<UInt> &);
AKANTU_GET_MACRO(MaterialList, material, const Vector<UInt> &);
protected:
Vector<UInt> material;
};
typedef FEMTemplate<IntegratorGauss,ShapeLagrange> MyFEMType;
SolidMechanicsModel(Mesh & mesh,
UInt spatial_dimension = 0,
const ID & id = "solid_mechanics_model",
const MemoryID & memory_id = 0);
virtual ~SolidMechanicsModel();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize completely the model
virtual void initFull(std::string material_file,
AnalysisMethod method = _explicit_dynamic);
/// initialize the fem object needed for boundary conditions
void initFEMBoundary(bool create_surface = true);
/// register the tags associated with the parallel synchronizer
void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL);
/// allocate all vectors
void initVectors();
/// initialize all internal arrays for materials
void initMaterials();
/// initialize the model
virtual void initModel();
/// init PBC synchronizer
void initPBC();
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/* ------------------------------------------------------------------------ */
/* PBC */
/* ------------------------------------------------------------------------ */
public:
/// change the equation number for proper assembly when using PBC
void changeEquationNumberforPBC(std::map<UInt,UInt> & pbc_pair);
/// synchronize Residual for output
void synchronizeResidual();
protected:
/// register PBC synchronizer
void registerPBCSynchronizer();
/* ------------------------------------------------------------------------ */
/* Explicit */
/* ------------------------------------------------------------------------ */
public:
/// initialize the stuff for the explicit scheme
void initExplicit();
/// initialize the array needed by updateResidual (residual, current_position)
void initializeUpdateResidualData();
/// update the current position vector
void updateCurrentPosition();
/// assemble the residual for the explicit scheme
void updateResidual(bool need_initialize = true);
/**
* \brief compute the acceleration from the residual
* this function is the explicit equivalent to solveDynamic in implicit
* In the case of lumped mass just divide the residual by the mass
* In the case of not lumped mass call solveDynamic<_acceleration_corrector>
*/
void updateAcceleration();
/// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix
void solveLumped(Vector<Real> & x,
const Vector<Real> & A,
const Vector<Real> & b,
const Vector<bool> & boundary,
Real alpha);
/// explicit integration predictor
void explicitPred();
/// explicit integration corrector
void explicitCorr();
/* ------------------------------------------------------------------------ */
/* Implicit */
/* ------------------------------------------------------------------------ */
public:
/// initialize the solver and the jacobian_matrix (called by initImplicit)
void initSolver(SolverOptions & options = _solver_no_options);
/// initialize the stuff for the implicit solver
void initImplicit(bool dynamic = false,
SolverOptions & solver_options = _solver_no_options);
/// solve Ma = f to get the initial acceleration
void initialAcceleration();
/// assemble the stiffness matrix
void assembleStiffnessMatrix();
/// solve @f[ A\delta u = f_ext - f_int @f] in displacement
void solveDynamic();
/// solve Ku = f
void solveStatic();
/// test the convergence (norm of increment)
bool testConvergenceIncrement(Real tolerance);
bool testConvergenceIncrement(Real tolerance, Real & error);
/// test the convergence (norm of residual)
bool testConvergenceResidual(Real tolerance);
bool testConvergenceResidual(Real tolerance, Real & error);
/// implicit time integration predictor
void implicitPred();
/// implicit time integration corrector
void implicitCorr();
protected:
/// finish the computation of residual to solve in increment
void updateResidualInternal();
/// compute the support reaction and store it in force
void updateSupportReaction();
/// compute A and solve @f[ A\delta u = f_ext - f_int @f]
template<NewmarkBeta::IntegrationSchemeCorrectorType type>
void solveDynamic(Vector<Real> & increment);
/* ------------------------------------------------------------------------ */
/* Explicit/Implicit */
/* ------------------------------------------------------------------------ */
public:
/// compute the stresses
void computeStresses();
/* ------------------------------------------------------------------------ */
/* Boundaries (solid_mechanics_model_boundary.cc) */
/* ------------------------------------------------------------------------ */
public:
class SurfaceLoadFunctor {
public:
virtual void traction(__attribute__ ((unused)) const types::Vector<Real> & position,
__attribute__ ((unused)) types::Vector<Real> & traction,
__attribute__ ((unused)) const types::Vector<Real> & normal,
__attribute__ ((unused)) Surface surface_id) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
virtual void stress(__attribute__ ((unused)) const types::Vector<Real> & position,
__attribute__ ((unused)) types::RMatrix & stress,
__attribute__ ((unused)) const types::Vector<Real> & normal,
__attribute__ ((unused)) Surface surface_id) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
};
/// compute force vector from a function(x,y,z) that describe stresses
void computeForcesFromFunction(BoundaryFunction in_function,
BoundaryFunctionType function_type) __attribute__((deprecated));
template<class Functor>
void computeForcesFromFunction(Functor & functor, BoundaryFunctionType function_type);
/// integrate a force on the boundary by providing a stress tensor
void computeForcesByStressTensor(const Vector<Real> & stresses,
const ElementType & type,
const GhostType & ghost_type);
/// integrate a force on the boundary by providing a traction vector
void computeForcesByTractionVector(const Vector<Real> & tractions,
const ElementType & type,
const GhostType & ghost_type);
/// synchronize the ghost element boundaries values
void synchronizeBoundaries();
/* ------------------------------------------------------------------------ */
/* Materials (solid_mechanics_model_material.cc) */
/* ------------------------------------------------------------------------ */
public:
/// register a material in the dynamic database
Material & registerNewMaterial(const ID & mat_type,
const std::string & opt_param = "");
template <typename M>
Material & registerNewCustomMaterial(const ID & mat_type,
const std::string & opt_param = "");
/// read the material files to instantiate all the materials
void readMaterials(const std::string & filename);
/// read a custom material with a keyword and class as template
template <typename M>
UInt readCustomMaterial(const std::string & filename,
const std::string & keyword);
/// Use a UIntData in the mesh to specify the material to use per element
void setMaterialIDsFromIntData(const std::string & data_name);
/* ------------------------------------------------------------------------ */
/* Mass (solid_mechanics_model_mass.cc) */
/* ------------------------------------------------------------------------ */
public:
/// assemble the lumped mass matrix
void assembleMassLumped();
/// assemble the mass matrix
void assembleMass();
protected:
/// assemble the lumped mass matrix for local and ghost elements
void assembleMassLumped(GhostType ghost_type);
void assembleMass(GhostType ghost_type);
/// fill a vector of rho
void computeRho(Vector<Real> & rho,
ElementType type,
GhostType ghost_type);
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
inline virtual UInt getNbDataForElements(const Vector<Element> & elements,
SynchronizationTag tag) const;
inline virtual void packElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag) const;
inline virtual void unpackElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag);
inline virtual UInt getNbDataToPack(SynchronizationTag tag) const;
inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const;
inline virtual void packData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag) const;
inline virtual void unpackData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag);
protected:
inline void splitElementByMaterial(const Vector<Element> & elements,
Vector<Element> * elements_per_mat) const;
inline virtual void packBarycenter(CommunicationBuffer & buffer,
const Vector<Element> & elements) const;
inline virtual void unpackBarycenter(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag);
/* ------------------------------------------------------------------------ */
/* Mesh Event Handler inherited members */
/* ------------------------------------------------------------------------ */
protected:
virtual void onNodesAdded (const Vector<UInt> & nodes_list,
const NewNodesEvent & event);
virtual void onNodesRemoved(const Vector<UInt> & element_list,
const Vector<UInt> & new_numbering,
const RemovedNodesEvent & event);
virtual void onElementsAdded (const Vector<Element> & nodes_list,
const NewElementsEvent & event);
virtual void onElementsRemoved(const Vector<Element> & element_list,
const ByElementTypeUInt & new_numbering,
const RemovedElementsEvent & event);
/* ------------------------------------------------------------------------ */
/* Dumpable interface */
/* ------------------------------------------------------------------------ */
public:
virtual void addDumpField(const std::string & field_id);
virtual void addDumpFieldVector(const std::string & field_id);
virtual void addDumpFieldTensor(const std::string & field_id);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// return the dimension of the system space
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
/// get the current value of the time step
AKANTU_GET_MACRO(TimeStep, time_step, Real);
/// set the value of the time step
AKANTU_SET_MACRO(TimeStep, time_step, Real);
/// get the value of the conversion from forces/ mass to acceleration
AKANTU_GET_MACRO(F_M2A, f_m2a, Real);
/// set the value of the conversion from forces/ mass to acceleration
AKANTU_SET_MACRO(F_M2A, f_m2a, Real);
/// get the SolidMechanicsModel::displacement vector
AKANTU_GET_MACRO(Displacement, *displacement, Vector<Real> &);
/// get the SolidMechanicsModel::current_position vector \warn only consistent
/// after a call to SolidMechanicsModel::updateCurrentPosition
AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector<Real> &);
/// get the SolidMechanicsModel::increment vector \warn only consistent if
/// SolidMechanicsModel::setIncrementFlagOn has been called before
AKANTU_GET_MACRO(Increment, *increment, Vector<Real> &);
/// get the lumped SolidMechanicsModel::mass vector
AKANTU_GET_MACRO(Mass, *mass, Vector<Real> &);
/// get the SolidMechanicsModel::velocity vector
AKANTU_GET_MACRO(Velocity, *velocity, Vector<Real> &);
/// get the SolidMechanicsModel::acceleration vector, updated by
/// SolidMechanicsModel::updateAcceleration
AKANTU_GET_MACRO(Acceleration, *acceleration, Vector<Real> &);
/// get the SolidMechanicsModel::force vector (boundary forces)
AKANTU_GET_MACRO(Force, *force, Vector<Real> &);
/// get the SolidMechanicsModel::residual vector, computed by
/// SolidMechanicsModel::updateResidual
AKANTU_GET_MACRO(Residual, *residual, Vector<Real> &);
/// get the SolidMechanicsModel::boundary vector
AKANTU_GET_MACRO(Boundary, *boundary, Vector<bool> &);
/// get the SolidMechnicsModel::incrementAcceleration vector
AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, Vector<Real> &);
/// get the value of the SolidMechanicsModel::increment_flag
AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool);
/// get the SolidMechanicsModel::element_material vector corresponding to a
/// given akantu::ElementType
// AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementMaterial, element_material, UInt);
// AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt);
//AKANTU_GET_MACRO(ElementMaterial, element_material, const ByElementTypeVector<UInt> &);
/// vectors containing local material element index for each global element index
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementIndexByMaterial, element_index_by_material, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementIndexByMaterial, element_index_by_material, UInt);
AKANTU_GET_MACRO(ElementIndexByMaterial, element_index_by_material, const ByElementTypeVector<UInt> &);
/// get a particular material
inline Material & getMaterial(UInt mat_index);
inline const Material & getMaterial(UInt mat_index) const;
/// give the number of materials
inline UInt getNbMaterials() const { return materials.size(); };
/// give the material internal index from its id
Int getInternalIndexFromID(const ID & id) const;
/// compute the stable time step
Real getStableTimeStep();
/// compute the potential energy
Real getPotentialEnergy();
/// compute the kinetic energy
Real getKineticEnergy();
/// compute the external work (for impose displacement, the velocity should be given too)
Real getExternalWork();
/// get the energies
Real getEnergy(const std::string & energy_id);
/// compute the energy for energy
Real getEnergy(const std::string & energy_id, ElementType & type, UInt index);
/// set the Contact object
AKANTU_SET_MACRO(Contact, contact, Contact *);
/**
* @brief set the SolidMechanicsModel::increment_flag to on, the activate the
* update of the SolidMechanicsModel::increment vector
*/
void setIncrementFlagOn();
/// get the stiffness matrix
AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &);
/// get the mass matrix
AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &);
inline FEM & getFEMBoundary(const ID & name = "");
/// get integrator
AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder &);
/// get access to the internal solver
AKANTU_GET_MACRO(Solver, *solver, Solver &);
protected:
/// compute the stable time step
Real getStableTimeStep(const GhostType & ghost_type);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// time step
Real time_step;
/// conversion coefficient form force/mass to acceleration
Real f_m2a;
/// displacements array
Vector<Real> * displacement;
/// lumped mass array
Vector<Real> * mass;
/// velocities array
Vector<Real> * velocity;
/// accelerations array
Vector<Real> * acceleration;
/// accelerations array
Vector<Real> * increment_acceleration;
/// forces array
Vector<Real> * force;
/// residuals array
Vector<Real> * residual;
/// boundaries array
Vector<bool> * boundary;
/// array of current position used during update residual
Vector<Real> * current_position;
/// mass matrix
SparseMatrix * mass_matrix;
/// velocity damping matrix
SparseMatrix * velocity_damping_matrix;
/// stiffness matrix
SparseMatrix * stiffness_matrix;
/// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta t^2}, d = \frac{\gamma}{\beta \Delta t} @f]
SparseMatrix * jacobian_matrix;
/// vectors containing local material element index for each global element index
ByElementTypeUInt element_index_by_material;
/// list of used materials
std::vector<Material *> materials;
/// integration scheme of second order used
IntegrationScheme2ndOrder * integrator;
/// increment of displacement
Vector<Real> * increment;
/// flag defining if the increment must be computed or not
bool increment_flag;
/// solver for implicit
Solver * solver;
/// object to resolve the contact
Contact * contact;
/// the spatial dimension
UInt spatial_dimension;
AnalysisMethod method;
Synchronizer * synch_parallel;
};
__END_AKANTU__
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "parser.hh"
#include "material.hh"
__BEGIN_AKANTU__
-
-#ifdef AKANTU_USE_IOHELPER
-# include "dumper_iohelper_tmpl_homogenizing_field.hh"
-# include "dumper_iohelper_tmpl_material_internal_field.hh"
-#endif
-
#include "solid_mechanics_model_tmpl.hh"
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "solid_mechanics_model_inline_impl.cc"
#endif
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */
diff --git a/src/synchronizer/grid_synchronizer.cc b/src/synchronizer/grid_synchronizer.cc
index 53f650fde..4b62c1167 100644
--- a/src/synchronizer/grid_synchronizer.cc
+++ b/src/synchronizer/grid_synchronizer.cc
@@ -1,480 +1,480 @@
/**
* @file grid_synchronizer.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Mon Oct 03 12:24:07 2011
*
* @brief implementation of the grid synchronizer
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "grid_synchronizer.hh"
#include "aka_grid_dynamic.hh"
#include "mesh.hh"
#include "fem.hh"
#include "static_communicator.hh"
#include "mesh_io.hh"
#include <iostream>
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
GridSynchronizer::GridSynchronizer(Mesh & mesh,
const ID & id,
MemoryID memory_id) :
DistributedSynchronizer(mesh, id, memory_id) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class E>
GridSynchronizer * GridSynchronizer::createGridSynchronizer(Mesh & mesh,
const SpatialGrid<E> & grid,
SynchronizerID id,
MemoryID memory_id) {
AKANTU_DEBUG_IN();
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
UInt nb_proc = comm.getNbProc();
UInt my_rank = comm.whoAmI();
GridSynchronizer & communicator = *(new GridSynchronizer(mesh, id, memory_id));
if(nb_proc == 1) return &communicator;
UInt spatial_dimension = mesh.getSpatialDimension();
Real * bounding_boxes = new Real[2 * spatial_dimension * nb_proc];
Real * my_bounding_box = bounding_boxes + 2 * spatial_dimension * my_rank;
// mesh.getLocalLowerBounds(my_bounding_box);
// mesh.getLocalUpperBounds(my_bounding_box + spatial_dimension);
const types::Vector<Real> & lower = grid.getLowerBounds();
const types::Vector<Real> & upper = grid.getUpperBounds();
for (UInt i = 0; i < spatial_dimension; ++i) {
my_bounding_box[i ] = lower(i);
my_bounding_box[spatial_dimension + i] = upper(i);
}
const types::Vector<Real> & spacing = grid.getSpacing();
AKANTU_DEBUG_INFO("Exchange of bounding box to detect the overlapping regions.");
comm.allGather(bounding_boxes, spatial_dimension * 2);
bool * intersects_proc = new bool[nb_proc];
std::fill_n(intersects_proc, nb_proc, true);
- UInt * first_cells = new UInt[3 * nb_proc];
- UInt * last_cells = new UInt[3 * nb_proc];
+ Int * first_cells = new Int[3 * nb_proc];
+ Int * last_cells = new Int[3 * nb_proc];
std::fill_n(first_cells, 3 * nb_proc, 0);
std::fill_n(first_cells, 3 * nb_proc, 0);
ByElementTypeUInt ** element_per_proc = new ByElementTypeUInt* [nb_proc];
for (UInt p = 0; p < nb_proc; ++p) element_per_proc[p] = NULL;
// check the overlapping between my box and the one from other processors
for (UInt p = 0; p < nb_proc; ++p) {
if(p == my_rank) continue;
Real * proc_bounding_box = bounding_boxes + 2 * spatial_dimension * p;
bool intersects = false;
- UInt * first_cell_p = first_cells + p * spatial_dimension;
- UInt * last_cell_p = last_cells + p * spatial_dimension;
+ Int * first_cell_p = first_cells + p * spatial_dimension;
+ Int * last_cell_p = last_cells + p * spatial_dimension;
for (UInt s = 0; s < spatial_dimension; ++s) {
// check overlapping of grid
intersects = Math::intersects(my_bounding_box[s],
my_bounding_box[spatial_dimension + s],
proc_bounding_box[s],
proc_bounding_box[spatial_dimension + s]);
intersects_proc[p] &= intersects;
if(intersects) {
AKANTU_DEBUG_INFO("I intersects with processor " << p << " in direction " << s);
// is point 1 of proc p in the dimension s in the range ?
bool point1 = Math::is_in_range(proc_bounding_box[s],
my_bounding_box[s],
my_bounding_box[s+spatial_dimension]);
// is point 2 of proc p in the dimension s in the range ?
bool point2 = Math::is_in_range(proc_bounding_box[s+spatial_dimension],
my_bounding_box[s],
my_bounding_box[s+spatial_dimension]);
Real start = 0.;
Real end = 0.;
if(point1 && !point2) {
/* |-----------| my_bounding_box(i)
* |-----------| proc_bounding_box(i)
* 1 2
*/
start = proc_bounding_box[s];
end = my_bounding_box[s+spatial_dimension];
} else if(point1 && point2) {
/* |-----------------| my_bounding_box(i)
* |-----------| proc_bounding_box(i)
* 1 2
*/
start = proc_bounding_box[s];
end = proc_bounding_box[s+spatial_dimension];
} else if(!point1 && point2) {
/* |-----------| my_bounding_box(i)
* |-----------| proc_bounding_box(i)
* 1 2
*/
start = my_bounding_box[s];
end = proc_bounding_box[s+spatial_dimension];
}
- first_cell_p[s] = grid.getCellID(start + spacing[s]/100., s);
- last_cell_p [s] = grid.getCellID(end - spacing[s]/100., s);
+ first_cell_p[s] = grid.getCellID(start - spacing[s]/100., s);
+ last_cell_p [s] = grid.getCellID(end + spacing[s]/100., s);
}
}
//create the list of cells in the overlapping
typedef typename SpatialGrid<E>::CellID CellID;
std::vector<CellID> * cell_ids = new std::vector<CellID>;
if(intersects_proc[p]) {
AKANTU_DEBUG_INFO("I intersects with processor " << p);
CellID cell_id(spatial_dimension);
// for (UInt i = 0; i < spatial_dimension; ++i) {
// if(first_cell_p[i] != 0) --first_cell_p[i];
// if(last_cell_p[i] != 0) ++last_cell_p[i];
// }
- for (UInt fd = first_cell_p[0]; fd <= last_cell_p[0]; ++fd) {
+ for (Int fd = first_cell_p[0]; fd <= last_cell_p[0]; ++fd) {
cell_id.setID(0, fd);
if(spatial_dimension == 1) {
cell_ids->push_back(cell_id);
}
else {
- for (UInt sd = first_cell_p[1]; sd <= last_cell_p[1] ; ++sd) {
+ for (Int sd = first_cell_p[1]; sd <= last_cell_p[1] ; ++sd) {
cell_id.setID(1, sd);
if(spatial_dimension == 2) {
cell_ids->push_back(cell_id);
} else {
- for (UInt ld = first_cell_p[2]; fd <= last_cell_p[2] ; ++ld) {
+ for (Int ld = first_cell_p[2]; fd <= last_cell_p[2] ; ++ld) {
cell_id.setID(2, ld);
cell_ids->push_back(cell_id);
}
}
}
}
}
// get the list of elements in the cells of the overlapping
typename std::vector<CellID>::iterator cur_cell_id = cell_ids->begin();
typename std::vector<CellID>::iterator last_cell_id = cell_ids->end();
std::set<Element> * to_send = new std::set<Element>();
for (; cur_cell_id != last_cell_id; ++cur_cell_id) {
typename SpatialGrid<E>::Cell::const_iterator cur_elem = grid.beginCell(*cur_cell_id);
typename SpatialGrid<E>::Cell::const_iterator last_elem = grid.endCell(*cur_cell_id);
for (; cur_elem != last_elem; ++cur_elem) {
to_send->insert(*cur_elem);
}
}
AKANTU_DEBUG_INFO("I have prepared " << to_send->size() << " elements to send to processor " << p);
delete cell_ids;
std::stringstream sstr; sstr << "element_per_proc_" << p;
element_per_proc[p] = new ByElementTypeUInt(sstr.str(), id);
ByElementTypeUInt & elempproc = *(element_per_proc[p]);
typename std::set<Element>::iterator elem = to_send->begin();
typename std::set<Element>::iterator last_elem = to_send->end();
for (; elem != last_elem; ++elem) {
ElementType type = elem->type;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
// /!\ this part must be slow due to the access in the ByElementTypeUInt
if(!elempproc.exists(type))
elempproc.alloc(0, nb_nodes_per_element, type, _not_ghost);
UInt global_connect[nb_nodes_per_element];
UInt * local_connect = mesh.getConnectivity(type).storage() + elem->element * nb_nodes_per_element;
for (UInt i = 0; i < nb_nodes_per_element; ++i) {
global_connect[i] = mesh.getNodeGlobalId(local_connect[i]);
}
elempproc(type).push_back(global_connect);
communicator.send_element[p].push_back(*elem);
}
delete to_send;
}
}
AKANTU_DEBUG_INFO("I have finished to compute intersection,"
<< " no it's time to communicate with my neighbors");
/**
* Sending loop, sends the connectivity asynchronously to all concerned proc
*/
std::vector<CommunicationRequest *> isend_requests;
for (UInt p = 0; p < nb_proc; ++p) {
if(p == my_rank) continue;
if(intersects_proc[p]) {
ByElementTypeUInt & elempproc = *(element_per_proc[p]);
ByElementTypeUInt::type_iterator it_type = elempproc.firstType(0, _not_ghost);
ByElementTypeUInt::type_iterator last_type = elempproc.lastType (0, _not_ghost);
UInt count = 0;
for (; it_type != last_type; ++it_type) {
Vector<UInt> & conn = elempproc(*it_type, _not_ghost);
UInt info[2];
info[0] = (UInt) *it_type;
info[1] = conn.getSize() * conn.getNbComponent();
AKANTU_DEBUG_INFO("I have " << conn.getSize() << " elements of type " << *it_type
<< " to send to processor " << p
<< " (communication tag : " << Tag::genTag(my_rank, count, DATA_TAG) << ")");
isend_requests.push_back(comm.asyncSend(info, 2, p, Tag::genTag(my_rank, count, SIZE_TAG)));
isend_requests.push_back(comm.asyncSend(conn.storage(),
conn.getSize() * conn.getNbComponent(),
p, Tag::genTag(my_rank, count, DATA_TAG)));
++count;
}
UInt info[2];
info[0] = (UInt) _not_defined;
info[1] = 0;
isend_requests.push_back(comm.asyncSend(info, 2, p, Tag::genTag(my_rank, count, SIZE_TAG)));
}
}
/**
* Receives the connectivity and store them in the ghosts elements
*/
Vector<UInt> & global_nodes_ids = const_cast<Vector<UInt> &>(mesh.getGlobalNodesIds());
Vector<Int> & nodes_type = const_cast<Vector<Int> &>(mesh.getNodesType());
std::vector<CommunicationRequest *> isend_nodes_requests;
UInt nb_nodes_to_recv[nb_proc];
UInt nb_total_nodes_to_recv = 0;
UInt nb_current_nodes = global_nodes_ids.getSize();
NewNodesEvent new_nodes;
NewElementsEvent new_elements;
Vector<UInt> * ask_nodes_per_proc = new Vector<UInt>[nb_proc];
for (UInt p = 0; p < nb_proc; ++p) {
nb_nodes_to_recv[p] = 0;
if(p == my_rank) continue;
Vector<UInt> & ask_nodes = ask_nodes_per_proc[p];
UInt count = 0;
if(intersects_proc[p]) {
ElementType type = _not_defined;
do {
UInt info[2] = { 0 };
comm.receive(info, 2, p, Tag::genTag(p, count, SIZE_TAG));
type = (ElementType) info[0];
if(type != _not_defined) {
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);;
UInt nb_element = info[1] / nb_nodes_per_element;
Vector<UInt> tmp_conn(nb_element, nb_nodes_per_element);
comm.receive(tmp_conn.storage(), info[1], p, Tag::genTag(p, count, DATA_TAG));
AKANTU_DEBUG_INFO("I will receive " << nb_element << " elements of type " << ElementType(info[0])
<< " from processor " << p
<< " (communication tag : " << Tag::genTag(p, count, DATA_TAG) << ")");
Vector<UInt> & ghost_connectivity = const_cast<Vector<UInt> &>(mesh.getConnectivity(type, _ghost));
UInt nb_ghost_element = ghost_connectivity.getSize();
Element element(type, 0, _ghost);
UInt conn[nb_nodes_per_element];
for (UInt el = 0; el < nb_element; ++el) {
UInt nb_node_to_ask_for_elem = 0;
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt gn = tmp_conn(el, n);
UInt ln = global_nodes_ids.find(gn);
if(ln == UInt(-1)) {
global_nodes_ids.push_back(gn);
nodes_type.push_back(-3); // pure ghost node
ln = nb_current_nodes;
new_nodes.getList().push_back(ln);
++nb_current_nodes;
ask_nodes.push_back(gn);
++nb_node_to_ask_for_elem;
}
conn[n] = ln;
}
// all the nodes are already known locally, the element should already exists
UInt c = UInt(-1);
if(nb_node_to_ask_for_elem == 0) {
c = ghost_connectivity.find(conn);
element.element = c;
}
if(c == UInt(-1)) {
element.element = nb_ghost_element;
++nb_ghost_element;
ghost_connectivity.push_back(conn);
new_elements.getList().push_back(element);
}
communicator.recv_element[p].push_back(element);
}
}
count++;
} while(type != _not_defined);
AKANTU_DEBUG_INFO("I have " << ask_nodes.getSize()
<< " missing nodes for elements coming from processor " << p
<< " (communication tag : " << Tag::genTag(my_rank, 0, ASK_NODES_TAG) << ")");
isend_nodes_requests.push_back(comm.asyncSend(ask_nodes.storage(), ask_nodes.getSize(),
p, Tag::genTag(my_rank, 0, ASK_NODES_TAG)));
nb_nodes_to_recv[p] = ask_nodes.getSize();
nb_total_nodes_to_recv += ask_nodes.getSize();
}
}
for (UInt p = 0; p < nb_proc; ++p) {
if(element_per_proc[p]) delete element_per_proc[p];
}
delete [] element_per_proc;
comm.waitAll(isend_requests);
comm.freeCommunicationRequest(isend_requests);
/**
* Sends requested nodes to proc
*/
Vector<Real> & nodes = const_cast<Vector<Real> &>(mesh.getNodes());
UInt nb_nodes = nodes.getSize();
std::vector<CommunicationRequest *> isend_coordinates_requests;
Vector<Real> * nodes_to_send_per_proc = new Vector<Real>[nb_proc];
for (UInt p = 0; p < nb_proc; ++p) {
if(p == my_rank || !intersects_proc[p]) continue;
Vector<UInt> asked_nodes;
CommunicationStatus status;
AKANTU_DEBUG_INFO("Waiting list of nodes to send to processor " << p
<< "(communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")");
comm.probe<UInt>(p, Tag::genTag(p, 0, ASK_NODES_TAG), status);
UInt nb_nodes_to_send = status.getSize();
asked_nodes.resize(nb_nodes_to_send);
AKANTU_DEBUG_INFO("I have " << nb_nodes_to_send
<< " nodes to send to processor " << p
<< " (communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")");
AKANTU_DEBUG_INFO("Getting list of nodes to send to processor " << p
<< " (communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")");
comm.receive(asked_nodes.storage(), nb_nodes_to_send, p, Tag::genTag(p, 0, ASK_NODES_TAG));
Vector<Real> & nodes_to_send = nodes_to_send_per_proc[p];
nodes_to_send.extendComponentsInterlaced(spatial_dimension, 1);
for (UInt n = 0; n < nb_nodes_to_send; ++n) {
UInt ln = global_nodes_ids.find(asked_nodes(n));
AKANTU_DEBUG_ASSERT(ln != UInt(-1), "The node [" << asked_nodes(n) << "] requested by proc " << p << " was not found locally!");
nodes_to_send.push_back(nodes.storage() + ln * spatial_dimension);
}
AKANTU_DEBUG_INFO("Sending the nodes to processor " << p
<< " (communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")");
isend_coordinates_requests.push_back(comm.asyncSend(nodes_to_send.storage(), nb_nodes_to_send * spatial_dimension, p, Tag::genTag(my_rank, 0, SEND_NODES_TAG)));
}
comm.waitAll(isend_nodes_requests);
comm.freeCommunicationRequest(isend_nodes_requests);
delete [] ask_nodes_per_proc;
// std::vector<CommunicationRequest *> irecv_nodes_requests;
nodes.resize(nb_total_nodes_to_recv + nb_nodes);
for (UInt p = 0; p < nb_proc; ++p) {
if((p != my_rank) && (nb_nodes_to_recv[p] > 0)) {
AKANTU_DEBUG_INFO("Receiving the nodes from processor " << p
<< " (communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")");
comm.receive(nodes.storage() + nb_nodes * spatial_dimension,
nb_nodes_to_recv[p] * spatial_dimension,
p, Tag::genTag(p, 0, SEND_NODES_TAG));
nb_nodes += nb_nodes_to_recv[p];
}
}
comm.waitAll(isend_coordinates_requests);
comm.freeCommunicationRequest(isend_coordinates_requests);
delete [] nodes_to_send_per_proc;
// comm.waitAll(irecv_nodes_requests);
// comm.freeCommunicationRequest(irecv_nodes_requests);
mesh.sendEvent(new_nodes);
mesh.sendEvent(new_elements);
AKANTU_DEBUG_OUT();
return &communicator;
}
/* -------------------------------------------------------------------------- */
template GridSynchronizer *
GridSynchronizer::createGridSynchronizer<QuadraturePoint>(Mesh & mesh,
const SpatialGrid<QuadraturePoint> & grid,
SynchronizerID id,
MemoryID memory_id);
template GridSynchronizer *
GridSynchronizer::createGridSynchronizer<Element>(Mesh & mesh,
const SpatialGrid<Element> & grid,
SynchronizerID id,
MemoryID memory_id);
__END_AKANTU__
diff --git a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d.cc b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d.cc
index c054c580c..f1976a4f5 100644
--- a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d.cc
+++ b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d.cc
@@ -1,175 +1,143 @@
/**
* @file test_heat_transfer_model_cube3d.cc
*
* @author Rui Wang <rui.wang@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
*
* @date Sun May 01 19:14:43 2011
*
* @brief test of the class HeatTransferModel on the 3d cube
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <fstream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "heat_transfer_model.hh"
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
#include "io_helper.hh"
#endif //AKANTU_USE_IOHELPER
using namespace akantu;
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
static void paraviewInit(HeatTransferModel & model, iohelper::Dumper & dumper);
static void paraviewDump(iohelper::Dumper & dumper);
iohelper::ElemType paraview_type = iohelper::TETRA1;
#endif
UInt spatial_dimension = 3;
ElementType type = _tetrahedron_4;
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[])
{
initialize(argc, argv);
Mesh mesh(spatial_dimension);
MeshIOMSH mesh_io;
mesh_io.read("cube1.msh", mesh);
HeatTransferModel model(mesh);
//initialize everything
model.initFull("material.dat");
//assemble the lumped capacity
model.assembleCapacityLumped();
//get and set stable time step
Real time_step = model.getStableTimeStep()*0.8;
std::cout << "time step is:" << time_step << std::endl;
model.setTimeStep(time_step);
/// boundary conditions
const Vector<Real> & nodes = mesh.getNodes();
Vector<bool> & boundary = model.getBoundary();
Vector<Real> & temperature = model.getTemperature();
UInt nb_nodes = mesh.getNbNodes();
//double t1, t2;
double length;
// t1 = 300.;
// t2 = 100.;
length = 1.;
for (UInt i = 0; i < nb_nodes; ++i) {
//temperature(i) = t1 - (t1 - t2) * sin(nodes(i, 0) * M_PI / length);
temperature(i) = 100.;
// if(nodes(i,0) < eps) {
// boundary(i) = true;
// temperature(i) = 100.0;
// }
//set the second boundary condition
// if(std::abs(nodes(i,0) - length) < eps) {
// boundary(i) = true;
// temperature(i) = 100.;
// }
//to insert a heat source
Real dx = nodes(i,0) - length/2.;
Real dy = nodes(i,1) - length/2.;
Real dz = nodes(i,2) - length/2.;
Real d = sqrt(dx*dx + dy*dy + dz*dz);
if(d < 0.1){
boundary(i) = true;
temperature(i) = 300.;
}
}
-#ifdef AKANTU_USE_IOHELPER
+
model.updateResidual();
- iohelper::DumperParaview dumper;
- paraviewInit(model,dumper);
-#endif
-
+ model.setBaseName("heat_transfer_cube3d");
+ model.addDumpField("temperature" );
+ model.addDumpField("temperature_rate");
+ model.addDumpField("residual" );
+ model.addDumpField("capacity_lumped" );
+ model.dump();
+
+
// //for testing
int max_steps = 1000;
for(int i=0; i<max_steps; i++)
{
model.explicitPred();
model.updateResidual();
model.explicitCorr();
-#ifdef AKANTU_USE_IOHELPER
- if(i % 100 == 0) paraviewDump(dumper);
-#endif
+
+ if(i % 100 == 0) model.dump();
if(i % 10 == 0) std::cout << "Step " << i << "/" << max_steps << std::endl;
}
std::cout << "Stable Time Step is : " << time_step << std::endl;
return 0;
}
-/* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_IOHELPER
-void paraviewInit(HeatTransferModel & model, iohelper::Dumper & dumper) {
- Mesh & mesh = model.getFEM().getMesh();
- UInt nb_nodes = mesh.getNbNodes();
- UInt nb_element = mesh.getNbElement(type);
-
-#pragma message "To change with new dumper"
- dumper.setMode(iohelper::TEXT);
- dumper.setPoints(mesh.getNodes().values,
- spatial_dimension, nb_nodes, "coordinates2");
- dumper.setConnectivity((int *) mesh.getConnectivity(type).values,
- paraview_type, nb_element, iohelper::C_MODE);
- // dumper.addNodeDataField(model.getTemperature().values,
- // 1, "temperature");
- // dumper.addNodeDataField(model.getTemperatureRate().values,
- // 1, "temperature_rate");
- // dumper.addNodeDataField(model.getResidual().values,
- // 1, "residual");
- // dumper.addNodeDataField(model.getCapacityLumped().values,
- // 1, "capacity_lumped");
- // dumper.addElemDataField(model.getTemperatureGradient(type).values,
- // spatial_dimension, "temperature_gradient");
- dumper.setPrefix("paraview/");
- dumper.init();
- dumper.dump();
-}
-
-/* -------------------------------------------------------------------------- */
-
-void paraviewDump(iohelper::Dumper & dumper) {
- dumper.dump();
-}
-#endif
-/* -------------------------------------------------------------------------- */
diff --git a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_istropic_conductivity.cc b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_istropic_conductivity.cc
index 7c53770eb..7aa820a95 100644
--- a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_istropic_conductivity.cc
+++ b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_istropic_conductivity.cc
@@ -1,150 +1,119 @@
/**
* @file test_heat_transfer_model_cube3d_istropic_conductivity.cc
*
* @author Rui Wang <rui.wang@epfl.ch>
*
* @date Sun May 01 19:14:43 2011
*
* @brief test of the class HeatTransferModel on the 3d cube
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "heat_transfer_model.hh"
#include <iostream>
#include <fstream>
#include <string.h>
using namespace std;
/* -------------------------------------------------------------------------- */
akantu::UInt spatial_dimension = 3;
-akantu:: ElementType type = akantu::_tetrahedron_4;
-/* -------------------------------------------------------------------------- */
-
-#ifdef AKANTU_USE_IOHELPER
-#include "io_helper.hh"
-iohelper::ElemType paraview_type = iohelper::TETRA1;
-
-static void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper);
-static void paraviewDump(iohelper::Dumper & dumper);
-#endif //AKANTU_USE_IOHELPER
-
-
-
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[])
{
akantu::initialize(argc, argv);
akantu::Mesh mesh(spatial_dimension);
akantu::MeshIOMSH mesh_io;
mesh_io.read("cube1.msh", mesh);
akantu::HeatTransferModel model(mesh);
//initialize everything
model.initFull("material.dat");
//assemble the lumped capacity
model.assembleCapacityLumped();
//get stable time step
akantu::Real time_step = model.getStableTimeStep()*0.8;
cout<<"time step is:"<<time_step<<endl;
model.setTimeStep(time_step);
/// boundary conditions
const akantu::Vector<akantu::Real> & nodes = model.getFEM().getMesh().getNodes();
akantu::Vector<bool> & boundary = model.getBoundary();
akantu::Vector<akantu::Real> & temperature = model.getTemperature();
akantu::Real eps = 1e-15;
double length = 1.;
akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
for (akantu::UInt i = 0; i < nb_nodes; ++i) {
//temperature(i) = t1 - (t1 - t2) * sin(nodes(i, 0) * M_PI / length);
temperature(i) = 100.;
if(nodes(i,0) < eps) {
boundary(i) = true;
temperature(i) = 300.;
}
//set the second boundary condition
if(std::abs(nodes(i,0) - length) < eps) {
boundary(i) = true;
temperature(i) = 300.;
}
// //to insert a heat source
// if(std::abs(nodes(i,0) - length/2.) < 0.025 && std::abs(nodes(i,1) - length/2.) < 0.025 && std::abs(nodes(i,2) - length/2.) < 0.025) {
// boundary(i) = true;
// temperature(i) = 300.;
// }
}
-#ifdef AKANTU_USE_IOHELPER
- iohelper::DumperParaview dumper;
- paraviewInit(model,dumper);
-#endif
+ model.updateResidual();
+ model.setBaseName("heat_transfer_cube3d_istropic_conductivity");
+ model.addDumpField("temperature" );
+ model.addDumpField("temperature_rate");
+ model.addDumpField("residual" );
+ model.addDumpField("capacity_lumped" );
+ model.dump();
// //for testing
int max_steps = 1000;
for(int i=0; i<max_steps; i++)
{
model.explicitPred();
model.updateResidual();
model.explicitCorr();
-#ifdef AKANTU_USE_IOHELPER
- if(i % 100 == 0)
- paraviewDump(dumper);
-#endif
+ if(i % 100 == 0) model.dump();
if(i % 10000 == 0)
std::cout << "Step " << i << "/" << max_steps << std::endl;
}
cout<< "\n\n Stable Time Step is : " << time_step << "\n \n" <<endl;
return 0;
}
-/* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_IOHELPER
-/* -------------------------------------------------------------------------- */
-
-void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper) {
-
-#pragma message "To change with new dumper"
-
-}
-/* -------------------------------------------------------------------------- */
-
-void paraviewDump(iohelper::Dumper & dumper) {
- // dumper.Dump();
-}
-/* -------------------------------------------------------------------------- */
-#endif
-/* -------------------------------------------------------------------------- */
-
diff --git a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_pbc.cc b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_pbc.cc
index ed5f633c1..1d4b55dae 100644
--- a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_pbc.cc
+++ b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_cube3d_pbc.cc
@@ -1,176 +1,137 @@
/**
* @file test_heat_transfer_model_cube3d_pbc.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Rui Wang <rui.wang@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
*
* @date Sun May 01 19:14:43 2011
*
* @brief test of the class HeatTransferModel on the 3d cube
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "heat_transfer_model.hh"
#include "pbc_synchronizer.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <fstream>
#include <string.h>
using namespace std;
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
#include "io_helper.hh"
#endif //AKANTU_USE_IOHELPER
#ifdef AKANTU_USE_IOHELPER
static void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper);
static void paraviewDump(iohelper::Dumper & dumper);
iohelper::ElemType paraview_type = iohelper::TETRA1;
#endif
akantu::UInt spatial_dimension = 3;
akantu:: ElementType type = akantu::_tetrahedron_4;
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[])
{
akantu::initialize(argc, argv);
//create mesh
akantu::Mesh mesh(spatial_dimension);
akantu::MeshIOMSH mesh_io;
mesh_io.read("cube_tet4.msh", mesh);
akantu::HeatTransferModel model(mesh);
//initialize everything
model.initFull("material.dat");
//initialize PBC
model.setPBC(1,1,1);
model.initPBC();
//assemble the lumped capacity
model.assembleCapacityLumped();
//get stable time step
akantu::Real time_step = model.getStableTimeStep()*0.8;
cout<<"time step is:"<<time_step<<endl;
model.setTimeStep(time_step);
// boundary conditions
const akantu::Vector<akantu::Real> & nodes = model.getFEM().getMesh().getNodes();
akantu::Vector<bool> & boundary = model.getBoundary();
akantu::Vector<akantu::Real> & temperature = model.getTemperature();
// double t1, t2;
double length;
// t1 = 300.;
// t2 = 100.;
length = 1.;
akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
for (akantu::UInt i = 0; i < nb_nodes; ++i) {
temperature(i) = 100.;
akantu::Real dx = nodes(i,0) - length/4.;
akantu::Real dy = nodes(i,1) - length/4.;
akantu::Real dz = nodes(i,2) - length/4.;
akantu::Real d = sqrt(dx*dx + dy*dy + dz*dz);
if(d < 0.1){
boundary(i) = true;
temperature(i) = 300.;
}
}
- /* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_IOHELPER
- iohelper::DumperParaview dumper;
- paraviewInit(model,dumper);
-#endif
+ model.updateResidual();
+ model.setBaseName("heat_transfer_cube3d_pbc");
+ model.addDumpField("temperature" );
+ model.addDumpField("temperature_rate");
+ model.addDumpField("residual" );
+ model.addDumpField("capacity_lumped" );
+ model.dump();
+
/* ------------------------------------------------------------------------ */
// //for testing
int max_steps = 1000;
/* ------------------------------------------------------------------------ */
for(int i=0; i<max_steps; i++)
{
model.explicitPred();
model.updateResidual();
model.explicitCorr();
-#ifdef AKANTU_USE_IOHELPER
- if(i % 100 == 0)
- paraviewDump(dumper);
-#endif
+ if(i % 100 == 0) model.dump();
+
if(i % 10 == 0)
std::cout << "Step " << i << "/" << max_steps << std::endl;
}
cout<< "\n\n Stable Time Step is : " << time_step << "\n \n" <<endl;
return 0;
}
-
-
-/* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_IOHELPER
-/* -------------------------------------------------------------------------- */
-
-void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper) {
- // akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
- // akantu::UInt nb_element = model.getFEM().getMesh().getNbElement(type);
-
-#pragma message "To change with new dumper"
- // dumper.SetMode(iohelper::TEXT);
- // dumper.SetPoints(model.getFEM().getMesh().getNodes().values,
- // spatial_dimension, nb_nodes, "coordinates2");
- // dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values,
- // paraview_type, nb_element, iohelper::C_MODE);
- // dumper.AddNodeDataField(model.getTemperature().values,
- // 1, "temperature");
- // dumper.AddNodeDataField(model.getResidual().values,
- // 1, "residual");
- // dumper.AddNodeDataField(model.getCapacityLumped().values,
- // 1, "capacity_lumped");
- // // dumper.AddElemDataField(model.getTemperatureGradient(type).values,
- // // spatial_dimension, "temperature_gradient");
- // // dumper.AddElemDataField(model.getbtkgt().values,
- // // 4, "btkgt");
- // dumper.SetPrefix("paraview/");
- // dumper.Init();
- // dumper.Dump();
-}
-
-/* -------------------------------------------------------------------------- */
-
-void paraviewDump(iohelper::Dumper & dumper) {
- // dumper.Dump();
-}
-
-/* -------------------------------------------------------------------------- */
-#endif
-/* -------------------------------------------------------------------------- */
diff --git a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc
index 498be9103..b31d0162d 100644
--- a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc
+++ b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc
@@ -1,169 +1,116 @@
/**
* @file test_heat_transfer_model_square2d.cc
*
*
* @date Sun May 01 19:14:43 2011
*
* @brief test of the class HeatTransferModel on the 3d cube
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "heat_transfer_model.hh"
#include "pbc_synchronizer.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <fstream>
#include <string>
using namespace std;
-/* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_IOHELPER
-# include "dumper_paraview.hh"
-#endif //AKANTU_USE_IOHELPER
-
-static void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper);
-static void paraviewDump(iohelper::Dumper & dumper);
-iohelper::ElemType paraview_type = iohelper::TRIANGLE1;
/* -------------------------------------------------------------------------- */
akantu::UInt spatial_dimension = 2;
-akantu:: ElementType type = akantu::_triangle_3;
-/* -------------------------------------------------------------------------- */
std::string base_name;
int main(int argc, char *argv[])
{
akantu::debug::setDebugLevel(akantu::dblWarning);
akantu::initialize(argc, argv);
//create mesh
akantu::Mesh mesh(spatial_dimension);
akantu::MeshIOMSH mesh_io;
mesh_io.read("square_tri3.msh", mesh);
-
+
akantu::HeatTransferModel model(mesh);
//initialize everything
model.initFull("material.dat");
-
+
//assemble the lumped capacity
model.assembleCapacityLumped();
-
+
//get stable time step
akantu::Real time_step = model.getStableTimeStep()*0.8;
cout<<"time step is:" << time_step << endl;
model.setTimeStep(time_step);
-
+
//boundary conditions
const akantu::Vector<akantu::Real> & nodes = model.getFEM().getMesh().getNodes();
akantu::Vector<bool> & boundary = model.getBoundary();
akantu::Vector<akantu::Real> & temperature = model.getTemperature();
double length;
length = 1.;
akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
for (akantu::UInt i = 0; i < nb_nodes; ++i) {
temperature(i) = 100.;
akantu::Real dx = nodes(i,0) - length/4.;
akantu::Real dy = 0.0;
akantu::Real dz = 0.0;
if (spatial_dimension > 1) dy = nodes(i,1) - length/4.;
if (spatial_dimension == 3) dz = nodes(i,2) - length/4.;
akantu::Real d = sqrt(dx*dx + dy*dy + dz*dz);
// if(dx < 0.0){
if(d < 0.1){
boundary(i) = true;
temperature(i) = 300.;
}
}
- iohelper::DumperParaview dumper;
- paraviewInit(model,dumper);
+ model.updateResidual();
+ model.setBaseName("heat_transfer_square2d");
+ model.addDumpField("temperature" );
+ model.addDumpField("temperature_rate");
+ model.addDumpField("residual" );
+ model.addDumpField("capacity_lumped" );
+ model.dump();
//main loop
- int max_steps = 1000;
+ int max_steps = 1500;
for(int i=0; i<max_steps; i++)
{
model.explicitPred();
model.updateResidual();
model.explicitCorr();
-
-#ifdef AKANTU_USE_IOHELPER
- if(i % 100 == 0)
- paraviewDump(dumper);
-#endif
+
+ if(i % 100 == 0) model.dump();
if(i % 10 == 0)
- std::cout << "Step " << i << "/" << max_steps << std::endl;
+ std::cout << "Step " << i << "/" << max_steps << std::endl;
}
cout<< "\n\n Stable Time Step is : " << time_step << "\n \n" <<endl;
return 0;
}
-
-void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper) {
- base_name = "coordinates2";
- akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
- akantu::UInt nb_element = model.getFEM().getMesh().getNbElement(type);
-
-#pragma message "To change with new dumper"
- // dumper.SetMode(iohelper::TEXT);
- dumper.setPoints(model.getFEM().getMesh().getNodes().values,
- spatial_dimension, nb_nodes, base_name);
- dumper.setConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values,
- paraview_type, nb_element, iohelper::C_MODE);
- dumper.addNodeDataField("temperature",model.getTemperature().values,
- 1,nb_nodes);
- // dumper.addNodeDataField(model.getTemperatureRate().values,
- // 1, "temperature_rate");
- // dumper.addNodeDataField(model.getResidual().values,
- // 1, "residual");
- // dumper.addNodeDataField(model.getCapacityLumped().values,
- // 1, "capacity_lumped");
- // dumper.addElemDataField(model.getTemperatureGradient(type).values,
- // spatial_dimension, "temperature_gradient");
-
- // dumper.addElemDataField(model.getConductivityOnQpoints(type).values,
- // spatial_dimension*spatial_dimension, "conductivity_qpoints");
-
-
- dumper.setPrefix("./");
- dumper.init();
-}
-
-/* -------------------------------------------------------------------------- */
-
-void paraviewDump(iohelper::Dumper & dumper) {
- static int dump_step = 0;
- std::stringstream sstr;
- sstr << base_name << "-";
- sstr.width(5);
- sstr.fill('0');
- sstr << dump_step;
-
- dumper.dump(sstr.str());
-
- ++dump_step;
-}
diff --git a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d_pbc.cc b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d_pbc.cc
index 43626810d..029d1e59d 100644
--- a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d_pbc.cc
+++ b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d_pbc.cc
@@ -1,172 +1,117 @@
/**
* @file test_heat_transfer_model_square2d_pbc.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Sun May 01 19:14:43 2011
*
* @brief test of the class HeatTransferModel on the 3d cube
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "heat_transfer_model.hh"
#include "pbc_synchronizer.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <fstream>
#include <string.h>
using namespace std;
/* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_IOHELPER
-#include "io_helper.hh"
-
-static void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper);
-static void paraviewDump(iohelper::Dumper & dumper);
-iohelper::ElemType paraview_type = iohelper::TRIANGLE1;
-
-#endif //AKANTU_USE_IOHELPER
-
-/* -------------------------------------------------------------------------- */
-akantu::UInt spatial_dimension = 2;
-akantu:: ElementType type = akantu::_triangle_3;
-/* -------------------------------------------------------------------------- */
-
int main(int argc, char *argv[])
{
+ akantu::UInt spatial_dimension = 2;
akantu::debug::setDebugLevel(akantu::dblWarning);
akantu::initialize(argc, argv);
//create mesh
akantu::Mesh mesh(spatial_dimension);
akantu::MeshIOMSH mesh_io;
mesh_io.read("square_tri3.msh", mesh);
-
+
akantu::HeatTransferModel model(mesh);
//initialize everything
model.initFull("material.dat");
-
+
//initialize PBC
model.setPBC(1,1,1);
model.initPBC();
-
+
//assemble the lumped capacity
model.assembleCapacityLumped();
-
+
//get stable time step
akantu::Real time_step = model.getStableTimeStep()*0.8;
cout<<"time step is:" << time_step << endl;
model.setTimeStep(time_step);
-
+
//boundary conditions
const akantu::Vector<akantu::Real> & nodes = model.getFEM().getMesh().getNodes();
akantu::Vector<bool> & boundary = model.getBoundary();
akantu::Vector<akantu::Real> & temperature = model.getTemperature();
double length;
length = 1.;
akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
for (akantu::UInt i = 0; i < nb_nodes; ++i) {
temperature(i) = 100.;
akantu::Real dx = nodes(i,0) - length/4.;
akantu::Real dy = 0.0;
akantu::Real dz = 0.0;
if (spatial_dimension > 1) dy = nodes(i,1) - length/4.;
if (spatial_dimension == 3) dz = nodes(i,2) - length/4.;
akantu::Real d = sqrt(dx*dx + dy*dy + dz*dz);
// if(dx < 0.0){
if(d < 0.1){
boundary(i) = true;
temperature(i) = 300.;
}
}
-#ifdef AKANTU_USE_IOHELPER
- iohelper::DumperParaview dumper;
- paraviewInit(model,dumper);
-#endif
+ model.updateResidual();
+ model.setBaseName("heat_transfer_square_2d_pbc");
+ model.addDumpField("temperature" );
+ model.addDumpField("temperature_rate");
+ model.addDumpField("residual" );
+ model.addDumpField("capacity_lumped" );
+ model.dump();
+
//main loop
int max_steps = 1000;
for(int i=0; i<max_steps; i++)
{
model.explicitPred();
model.updateResidual();
model.explicitCorr();
-
-#ifdef AKANTU_USE_IOHELPER
- if(i % 100 == 0)
- paraviewDump(dumper);
-#endif
+ if(i % 100 == 0) model.dump();
if(i % 10 == 0)
std::cout << "Step " << i << "/" << max_steps << std::endl;
}
cout<< "\n\n Stable Time Step is : " << time_step << "\n \n" <<endl;
return 0;
}
-/* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_IOHELPER
-/* -------------------------------------------------------------------------- */
-
-void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper) {
- // akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
- // akantu::UInt nb_element = model.getFEM().getMesh().getNbElement(type);
-
-#pragma message "To change with new dumper"
- // // dumper.SetMode(iohelper::TEXT);
- // dumper.SetPoints(model.getFEM().getMesh().getNodes().values,
- // spatial_dimension, nb_nodes, "coordinates2");
- // dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values,
- // paraview_type, nb_element, iohelper::C_MODE);
- // dumper.AddNodeDataField(model.getTemperature().values,
- // 1, "temperature");
- // dumper.AddNodeDataField(model.getTemperatureRate().values,
- // 1, "temperature_rate");
- // dumper.AddNodeDataField(model.getResidual().values,
- // 1, "residual");
- // dumper.AddNodeDataField(model.getCapacityLumped().values,
- // 1, "capacity_lumped");
- // dumper.AddElemDataField(model.getTemperatureGradient(type).values,
- // spatial_dimension, "temperature_gradient");
-
- // dumper.AddElemDataField(model.getConductivityOnQpoints(type).values,
- // spatial_dimension*spatial_dimension, "conductivity_qpoints");
-
-
- // dumper.SetPrefix("paraview/");
- // dumper.Init();
-}
-
-/* -------------------------------------------------------------------------- */
-
-void paraviewDump(iohelper::Dumper & dumper) {
- // dumper.Dump();
-}
-
-/* -------------------------------------------------------------------------- */
-#endif
-/* -------------------------------------------------------------------------- */
diff --git a/test/test_synchronizer/test_grid_synchronizer.cc b/test/test_synchronizer/test_grid_synchronizer.cc
index 8bb041b20..f92559127 100644
--- a/test/test_synchronizer/test_grid_synchronizer.cc
+++ b/test/test_synchronizer/test_grid_synchronizer.cc
@@ -1,331 +1,334 @@
/**
* @file test_grid_synchronizer.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Fri Nov 25 17:00:17 2011
*
* @brief test the GridSynchronizer object
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_grid_dynamic.hh"
#include "mesh.hh"
-#include "mesh_io.hh"
#include "grid_synchronizer.hh"
#include "mesh_partition.hh"
#include "synchronizer_registry.hh"
#ifdef AKANTU_USE_IOHELPER
# include "io_helper.hh"
#endif //AKANTU_USE_IOHELPER
using namespace akantu;
+typedef std::map<std::pair<Element, Element>, Real> pair_list;
+
+static void updatePairList(const ByElementTypeReal & barycenter,
+ const SpatialGrid<Element> & grid,
+ Real radius,
+ pair_list & neighbors) {
+ AKANTU_DEBUG_IN();
+
+ GhostType ghost_type = _not_ghost;
+
+ Element e;
+ e.ghost_type = ghost_type;
+
+ // generate the pair of neighbor depending of the cell_list
+ ByElementTypeReal::type_iterator it = barycenter.firstType(0, ghost_type);
+ ByElementTypeReal::type_iterator last_type = barycenter.lastType(0, ghost_type);
+ for(; it != last_type; ++it) {
+ // loop over quad points
+
+ e.type = *it;
+ e.element = 0;
+
+ const Vector<Real> & barycenter_vect = barycenter(*it, ghost_type);
+ UInt sp = barycenter_vect.getNbComponent();
+
+ Vector<Real>::const_iterator< types::Vector<Real> > bary =
+ barycenter_vect.begin(sp);
+ Vector<Real>::const_iterator< types::Vector<Real> > bary_end =
+ barycenter_vect.end(sp);
+
+ for(;bary != bary_end; ++bary, e.element++) {
+ SpatialGrid<Element>::CellID cell_id = grid.getCellID(*bary);
+ SpatialGrid<Element>::neighbor_cells_iterator first_neigh_cell =
+ grid.beginNeighborCells(cell_id);
+ SpatialGrid<Element>::neighbor_cells_iterator last_neigh_cell =
+ grid.endNeighborCells(cell_id);
+
+ // loop over neighbors cells of the one containing the current element
+ for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) {
+ SpatialGrid<Element>::Cell::const_iterator first_neigh_el =
+ grid.beginCell(*first_neigh_cell);
+ SpatialGrid<Element>::Cell::const_iterator last_neigh_el =
+ grid.endCell(*first_neigh_cell);
+
+ // loop over the quadrature point in the current cell of the cell list
+ for (;first_neigh_el != last_neigh_el; ++first_neigh_el){
+ const Element & elem = *first_neigh_el;
+
+ Vector<Real>::const_iterator< types::Vector<Real> > neigh_it =
+ barycenter(elem.type, elem.ghost_type).begin(sp);
+
+ const types::RVector & neigh_bary = neigh_it[elem.element];
+
+ Real distance = bary->distance(neigh_bary);
+ if(distance <= radius) {
+ std::pair<Element, Element> pair = std::make_pair(e, elem);
+ pair_list::iterator p = neighbors.find(pair);
+ if(p != neighbors.end()) {
+ AKANTU_DEBUG_ERROR("Pair already registered [" << e << " " << elem << "] -> " << p->second << " " << distance);
+ } else {
+ neighbors[pair] = distance;
+ }
+ }
+ }
+ }
+ }
+ }
+ AKANTU_DEBUG_OUT();
+}
+
class TestAccessor : public DataAccessor {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
- TestAccessor(const Mesh & mesh);
- ~TestAccessor();
+ TestAccessor(const Mesh & mesh, const ByElementTypeReal & barycenters);
- AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Barycenter, barycenter, Real);
+ AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Barycenter, barycenters, Real);
/* ------------------------------------------------------------------------ */
/* Ghost Synchronizer inherited members */
/* ------------------------------------------------------------------------ */
protected:
virtual UInt getNbDataForElements(const Vector<Element> & elements,
SynchronizationTag tag) const;
virtual void packElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag) const;
virtual void unpackElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
- std::string id;
-
- ByElementTypeReal barycenter;
-
+ const ByElementTypeReal & barycenters;
const Mesh & mesh;
};
/* -------------------------------------------------------------------------- */
/* TestSynchronizer implementation */
/* -------------------------------------------------------------------------- */
-TestAccessor::TestAccessor(const Mesh & mesh) : barycenter("barycenter", id), mesh(mesh) {
- UInt spatial_dimension = mesh.getSpatialDimension();
-
- id = "test_synchronizer";
-
- Mesh::ConnectivityTypeList::const_iterator it;
- try {
- const Mesh::ConnectivityTypeList & ghost_type_list = mesh.getConnectivityTypeList(_ghost);
- for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) {
- if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
-
- UInt nb_ghost_element = mesh.getNbElement(*it,_ghost);
- barycenter.alloc(nb_ghost_element, spatial_dimension, *it);
- }
- } catch (debug::Exception & e) { std::cout << e << std::endl; }
-}
-
-TestAccessor::~TestAccessor() {
-
-}
+TestAccessor::TestAccessor(const Mesh & mesh,
+ const ByElementTypeReal & barycenters) : barycenters(barycenters), mesh(mesh) { }
UInt TestAccessor::getNbDataForElements(const Vector<Element> & elements,
__attribute__ ((unused)) SynchronizationTag tag) const {
- return Mesh::getSpatialDimension(elements(0).type) * sizeof(Real) * elements.getSize();
+ if(elements.getSize())
+ return Mesh::getSpatialDimension(elements(0).type) * sizeof(Real) * elements.getSize();
+ else
+ return 0;
}
void TestAccessor::packElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
__attribute__ ((unused)) SynchronizationTag tag) const {
UInt spatial_dimension = mesh.getSpatialDimension();
Vector<Element>::const_iterator<Element> bit = elements.begin();
Vector<Element>::const_iterator<Element> bend = elements.end();
for (; bit != bend; ++bit) {
const Element & element = *bit;
- types::RVector bary(barycenter(element.type,
- element.ghost_type).storage() + element.element * spatial_dimension,
- spatial_dimension);;
- buffer << barycenter;
+ types::RVector bary(this->barycenters(element.type, element.ghost_type).storage()
+ + element.element * spatial_dimension,
+ spatial_dimension);
+ buffer << bary;
}
}
void TestAccessor::unpackElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
__attribute__ ((unused)) SynchronizationTag tag) {
UInt spatial_dimension = mesh.getSpatialDimension();
Vector<Element>::const_iterator<Element> bit = elements.begin();
Vector<Element>::const_iterator<Element> bend = elements.end();
for (; bit != bend; ++bit) {
const Element & element = *bit;
- types::RVector barycenter_loc(spatial_dimension);
- mesh.getBarycenter(element.element, element.type, barycenter_loc.storage(), element.ghost_type);
+ types::RVector barycenter_loc(this->barycenters(element.type,element.ghost_type).storage()
+ + element.element * spatial_dimension,
+ spatial_dimension);
- types::RVector barycenter(spatial_dimension);
- buffer >> barycenter;
+ types::RVector bary(spatial_dimension);
+ buffer >> bary;
Real tolerance = 1e-15;
for (UInt i = 0; i < spatial_dimension; ++i) {
- if(!(std::abs(barycenter(i) - barycenter_loc(i)) <= tolerance))
+ if(!(std::abs(bary(i) - barycenter_loc(i)) <= tolerance))
AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: "
<< element
<< "(barycenter[" << i << "] = " << barycenter_loc(i)
- << " and buffer[" << i << "] = " << barycenter(i) << ") - tag: " << tag);
+ << " and buffer[" << i << "] = " << bary(i) << ") - tag: " << tag);
}
}
}
/* -------------------------------------------------------------------------- */
/* Main */
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[]) {
akantu::initialize(argc, argv);
UInt spatial_dimension = 2;
- ElementType type = _triangle_3;
+ Real radius = 0.1;
Mesh mesh(spatial_dimension);
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
Int psize = comm.getNbProc();
Int prank = comm.whoAmI();
if(prank == 0) {
- MeshIOMSH mesh_io;
- mesh_io.read("triangle.msh", mesh);
+ mesh.read("triangle.msh");
MeshPartition * partition = new MeshPartitionScotch(mesh, spatial_dimension);
partition->partitionate(psize);
DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition);
delete partition;
} else {
DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL);
}
-#ifdef AKANTU_USE_IOHELPER
- double * part;
- unsigned int nb_nodes = mesh.getNbNodes();
- unsigned int nb_element = mesh.getNbElement(type);
-
- iohelper::DumperParaview dumper;
- dumper.setMode(iohelper::TEXT);
- dumper.setParallelContext(prank, psize);
- dumper.setPoints(mesh.getNodes().values, spatial_dimension, nb_nodes, "test-grid-synchronizer");
- dumper.setConnectivity((int*) mesh.getConnectivity(type).values,
- iohelper::TRIANGLE1, nb_element, iohelper::C_MODE);
- part = new double[nb_element];
- for (unsigned int i = 0; i < nb_element; ++i)
- part[i] = prank;
- dumper.addElemDataField("partitions", part, 1, nb_element);
- dumper.setPrefix("paraview/");
- dumper.init();
- dumper.dump();
- delete [] part;
-
- iohelper::DumperParaview dumper_ghost;
- dumper_ghost.setMode(iohelper::TEXT);
- dumper_ghost.setParallelContext(prank, psize);
-
- try {
- unsigned int nb_ghost_element = mesh.getNbElement(type,_ghost);
- dumper_ghost.setPoints(mesh.getNodes().values, spatial_dimension, nb_nodes, "test-grid-synchronizer-ghost");
- dumper_ghost.setConnectivity((int*) mesh.getConnectivity(type,_ghost).values,
- iohelper::TRIANGLE1, nb_ghost_element, iohelper::C_MODE);
- part = new double[nb_ghost_element];
- for (unsigned int i = 0; i < nb_ghost_element; ++i)
- part[i] = prank;
- dumper_ghost.addElemDataField("partitions", part, 1, nb_ghost_element);
- dumper_ghost.setPrefix("paraview/");
- dumper_ghost.init();
- dumper_ghost.dump();
- delete [] part;
- } catch (debug::Exception & e) { std::cout << e << std::endl; }
-#endif
-
-
- comm.barrier();
-
mesh.computeBoundingBox();
Real lower_bounds[spatial_dimension];
Real upper_bounds[spatial_dimension];
mesh.getLowerBounds(lower_bounds);
mesh.getUpperBounds(upper_bounds);
types::Vector<Real> spacing(spatial_dimension);
types::Vector<Real> center(spatial_dimension);
for (UInt i = 0; i < spatial_dimension; ++i) {
- spacing[i] = (upper_bounds[i] - lower_bounds[i]) / 10.;
+ spacing[i] = radius * 1.2;
center[i] = (upper_bounds[i] + lower_bounds[i]) / 2.;
}
- SpacingGrid<Element> grid(spatial_dimension, spacing, center);
+ SpatialGrid<Element> grid(spatial_dimension, spacing, center);
GhostType ghost_type = _not_ghost;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
ByElementTypeReal barycenters("", "", 0);
- mesh.initByElementTypeVector(barycenters, spatial_dimension, 0);
+ mesh.initByElementTypeVector(barycenters, spatial_dimension, spatial_dimension);
Element e;
e.ghost_type = ghost_type;
- it = mesh.firstType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
- UInt nb_element = mesh.getNbElement(*it);
+ UInt nb_element = mesh.getNbElement(*it, ghost_type);
e.type = *it;
- Vector<Real> & barycenter = barycenters(*it);
+ Vector<Real> & barycenter = barycenters(*it, ghost_type);
barycenter.resize(nb_element);
Vector<Real>::iterator<types::RVector> bary_it = barycenter.begin(spatial_dimension);
for (UInt elem = 0; elem < nb_element; ++elem) {
- mesh.getBarycenter(elem, *it, bary_it->storage());
+ mesh.getBarycenter(elem, *it, bary_it->storage(), ghost_type);
e.element = elem;
grid.insert(e, *bary_it);
++bary_it;
}
}
- MeshIOMSH mesh_io;
-
std::stringstream sstr; sstr << "mesh_" << prank << ".msh";
- mesh_io.write(sstr.str(), mesh);
+ mesh.write(sstr.str());
+
+ std::cout << "Pouet 1" << std::endl;
+
+ GridSynchronizer * grid_communicator = GridSynchronizer::createGridSynchronizer(mesh, grid);
+
+ std::cout << "Pouet 2" << std::endl;
+
+ ghost_type = _ghost;
+
+ it = mesh.firstType(spatial_dimension, ghost_type);
+ last_type = mesh.lastType(spatial_dimension, ghost_type);
+ e.ghost_type = ghost_type;
+ for(; it != last_type; ++it) {
+ UInt nb_element = mesh.getNbElement(*it, ghost_type);
+ e.type = *it;
+ Vector<Real> & barycenter = barycenters(*it, ghost_type);
+ barycenter.resize(nb_element);
+
+ Vector<Real>::iterator<types::RVector> bary_it = barycenter.begin(spatial_dimension);
+ for (UInt elem = 0; elem < nb_element; ++elem) {
+ mesh.getBarycenter(elem, *it, bary_it->storage(), ghost_type);
+ e.element = elem;
+ grid.insert(e, *bary_it);
+ ++bary_it;
+ }
+ }
Mesh grid_mesh(spatial_dimension, "grid_mesh", 0);
std::stringstream sstr_grid; sstr_grid << "grid_mesh_" << prank << ".msh";
grid.saveAsMesh(grid_mesh);
- mesh_io.write(sstr_grid.str(), grid_mesh);
+ grid_mesh.write(sstr_grid.str());
- GridSynchronizer * grid_communicator = GridSynchronizer::createGridSynchronizer(mesh, grid);
+ std::cout << "Pouet 3" << std::endl;
+
+ pair_list neighbors;
+ updatePairList(barycenters, grid, radius, neighbors);
+ pair_list::iterator nit = neighbors.begin();
+ pair_list::iterator nend = neighbors.end();
+
+ std::stringstream sstrp; sstrp << "pairs_" << prank;
+ std::ofstream fout(sstrp.str().c_str());
+ for(;nit != nend; ++nit) {
+ fout << "[" << nit->first.first << "," << nit->first.second << "] -> "
+ << nit->second << std::endl;
+ }
+
+ fout.close();
+
+ std::cout << "Pouet 4" << std::endl;
AKANTU_DEBUG_INFO("Creating TestAccessor");
- TestAccessor test_accessor(mesh);
+ TestAccessor test_accessor(mesh, barycenters);
SynchronizerRegistry synch_registry(test_accessor);
synch_registry.registerSynchronizer(*grid_communicator, _gst_test);
AKANTU_DEBUG_INFO("Synchronizing tag");
synch_registry.synchronize(_gst_test);
-
-#ifdef AKANTU_USE_IOHELPER
- try {
- UInt nb_ghost_element = mesh.getNbElement(type, _ghost);
- UInt nb_nodes = mesh.getNbNodes();
-
- iohelper::DumperParaview dumper_ghost;
- dumper_ghost.setMode(iohelper::TEXT);
- dumper_ghost.setParallelContext(prank, psize);
-
- dumper_ghost.setMode(iohelper::TEXT);
- dumper_ghost.setParallelContext(prank, psize);
- dumper_ghost.setPoints(mesh.getNodes().values, spatial_dimension, nb_nodes, "test-grid-synchronizer-ghost");
- dumper_ghost.setConnectivity((int*) mesh.getConnectivity(type,_ghost).values,
- iohelper::TRIANGLE1, nb_ghost_element, iohelper::C_MODE);
- part = new double[nb_ghost_element];
- for (unsigned int i = 0; i < nb_ghost_element; ++i)
- part[i] = prank;
- dumper_ghost.addElemDataField("partitions", part, 1, nb_ghost_element);
- dumper_ghost.dump();
- delete [] part;
-
- unsigned int nb_grid_element = grid_mesh.getNbElement(_quadrangle_4);
- iohelper::DumperParaview dumper_grid;
- dumper_grid.setMode(iohelper::TEXT);
- dumper_grid.setParallelContext(prank, psize);
- dumper_grid.setPoints(grid_mesh.getNodes().values, spatial_dimension,
- grid_mesh.getNbNodes(), "test-grid-synchronizer-grid");
- dumper_grid.setConnectivity((int*) grid_mesh.getConnectivity(_quadrangle_4).values,
- iohelper::QUAD1, nb_grid_element, iohelper::C_MODE);
-
- part = new double[nb_grid_element];
- std::fill_n(part, nb_grid_element, prank);
-
- dumper_grid.addElemDataField("partitions", part, 1, nb_grid_element);
- dumper_grid.setPrefix("paraview/");
- dumper_grid.init();
- dumper_grid.dump();
- delete [] part;
- } catch(debug::Exception & e) { std::cout << e << std::endl; }
-#endif //AKANTU_USE_IOHELPER
-
-
akantu::finalize();
return EXIT_SUCCESS;
}