diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh index b993cfc3c..fb260a805 100644 --- a/src/common/aka_types.hh +++ b/src/common/aka_types.hh @@ -1,509 +1,509 @@ /** * @file aka_types.hh * @author Nicolas Richart * @date Wed Feb 16 20:28:13 2011 * * @brief description of the "simple" 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_error.hh" #include "aka_math.hh" #include "aka_vector.hh" /* -------------------------------------------------------------------------- */ #include #ifndef __INTEL_COMPILER #include #else #include #endif /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_TYPES_HH__ #define __AKANTU_AKA_TYPES_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* maps */ /* -------------------------------------------------------------------------- */ #ifndef __INTEL_COMPILER template struct unordered_map { typedef typename std::tr1::unordered_map type; }; #else template struct unordered_map { typedef typename std::map type; }; #endif namespace types { class Matrix; /* ------------------------------------------------------------------------ */ /* Vector */ /* ------------------------------------------------------------------------ */ template class Vector { public: Vector() : n(0), values(NULL), wrapped(false) {} Vector(UInt n, T def = T()) : n(n), values(new T[n]), wrapped(false) { std::fill_n(values, n, def); } Vector(T* data, UInt n) : n(n), values(data), wrapped(true) {} Vector(const Vector & src) { wrapped = src.wrapped; n = src.n; if (src.wrapped) { values = src.values; } else { values = new T[n]; memcpy(this->values, src.values, n * sizeof(T)); } } virtual ~Vector() { if(!wrapped) delete [] values; }; /* ---------------------------------------------------------------------- */ UInt size() const { return n; } T * storage() const { return values; } /* ---------------------------------------------------------------------- */ void shallowCopy(const Vector & src) { if(!wrapped) delete [] values; this->n = src.n; this->wrapped = true; this->values = src.values; } /* ---------------------------------------------------------------------- */ inline T& operator()(UInt i) { return *(values + i); }; inline const T& operator()(UInt i) const { return *(values + i); }; inline T& operator[](UInt idx) { return *(values + idx); }; inline const T& operator[](UInt idx) const { return *(values + idx); }; /* ---------------------------------------------------------------------- */ inline Vector & operator=(const Vector & src) { if(this != &src) { if (wrapped) { AKANTU_DEBUG_ASSERT(n == src.n, "vectors of different size"); memcpy(this->values, src.values, n * sizeof(T)); } else { n = src.n; delete []values; values = new T[n]; memcpy(this->values, src.values, n * sizeof(T)); } } return *this; } /* ---------------------------------------------------------------------- */ inline Vector & operator+=(const Vector & vect) { T * a = this->storage(); T * b = vect.storage(); for (UInt i = 0; i < n; ++i) *(a++) += *(b++); return *this; } /* ---------------------------------------------------------------------- */ inline Vector & operator+=(const T & x) { T * a = this->values; for (UInt i = 0; i < n; ++i) *(a++) += x; return *this; } /* ---------------------------------------------------------------------- */ inline Vector & operator-=(const Vector & vect) { T * a = this->storage(); T * b = vect.storage(); for (UInt i = 0; i < n; ++i) *(a++) -= *(b++); return *this; } /* ---------------------------------------------------------------------- */ inline Vector & operator*=(const T & scalar) { T * a = this->storage(); for (UInt i = 0; i < n; ++i) *(a++) *= scalar; return *this; } /* -------------------------------------------------------------------------- */ inline Vector & operator=(T & scalar) { T * a = this->storage(); for (UInt i = 0; i < n; ++i) *(a++) = scalar; return *this; } /* ---------------------------------------------------------------------- */ inline Vector & operator/=(const T & x) { T * a = this->values; for (UInt i = 0; i < n; ++i) *(a++) /= x; return *this; } /* ---------------------------------------------------------------------- */ inline Real dot(const Vector & vect) { return Math::vectorDot(values, vect.storage(), n); } /* ---------------------------------------------------------------------- */ inline Vector & crossProduct(const Vector & v1, const Vector & v2) { AKANTU_DEBUG_ASSERT(n == 3, "crossProduct is only defined in 3D"); AKANTU_DEBUG_ASSERT(n == v1.n && n == v2.n, "crossProduct is not a valid operation non matching size vectors"); // for (UInt i = 0; i < n; ++i) { // values[i] = // v1((i+1) % n) * v2((i+2) % n) - // v1((i+2) % n) * v1((i+1) % n); // } Math::vectorProduct3(v1.values, v2.values, this->values); return *this; } /* ---------------------------------------------------------------------- */ inline void clear() { memset(values, 0, n * sizeof(T)); }; template inline void mul(const Matrix & A, const Vector & x, Real alpha = 1.0); /* ---------------------------------------------------------------------- */ inline Real norm() const { return Math::norm(this->n, this->values); } /* ---------------------------------------------------------------------- */ /// norm of (*this - x) inline Real distance(const Vector & y) const { Real * vx = values; Real * vy = y.storage(); Real sum_2 = 0; for (UInt i = 0; i < n; ++i, ++vx, ++vy) sum_2 += (*vx - *vy)*(*vx - *vy); return sqrt(sum_2); } /* ---------------------------------------------------------------------- */ /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "types::Vector<" << debug::demangle(typeid(T).name()) << "> [" << n <<"] :" << std::endl; stream << space << AKANTU_INDENT << "| "; for (UInt i = 0; i < n; ++i) { stream << values[i] << " "; } stream << "|" << std::endl; } friend class ::akantu::Vector; protected: UInt n; T * values; bool wrapped; }; typedef Vector RVector; // support operations for the creation of other vectors template Vector operator*(T scalar, const Vector& a); template Vector operator+(const Vector& a, const Vector& b); template Vector operator-(const Vector& a, const Vector& b); /* -------------------------------------------------------------------------- */ template Vector operator*(T scalar, const Vector& a) { Vector r = a; r *= scalar; return r; } template Vector operator+(const Vector& a, const Vector& b) { Vector r = a; r += b; return r; } template Vector operator-(const Vector& a, const Vector& b) { Vector r = a; r -= b; return r; } /* ------------------------------------------------------------------------ */ /* Matrix */ /* ------------------------------------------------------------------------ */ class Matrix { public: Matrix() : m(0), n(0), values(NULL), wrapped(true) {}; Matrix(UInt m, UInt n, Real def = 0) : m(m), n(n), values(new Real[n*m]), wrapped(false) { std::fill_n(values, n*m, def); }; Matrix(Real * data, UInt m, UInt n) : m(m), n(n), values(data), wrapped(true) {}; Matrix(const Matrix & src) { m = src.m; n = src.n; values = src.values; wrapped = src.wrapped; const_cast(src).wrapped = true; }; virtual ~Matrix() { if(!wrapped) delete [] values; }; /* ---------------------------------------------------------------------- */ UInt size() const { return n*m; }; UInt rows() const { return m; }; UInt cols() const { return n; }; Real * storage() const { return values; }; /* ---------------------------------------------------------------------- */ void shallowCopy(const Matrix & src) { if(!wrapped) delete [] values; this->n = src.n; this->m = src.m; this->wrapped = true; this->values = src.values; } /* ---------------------------------------------------------------------- */ inline Real& operator()(UInt i, UInt j) { return *(values + i*n + j); }; inline const Real& operator()(UInt i, UInt j) const { return *(values + i*n + j); }; inline Real& operator[](UInt idx) { return *(values + idx); }; inline const Real& operator[](UInt idx) const { return *(values + idx); }; inline Matrix & operator=(const Matrix & mat) { if(this != &mat) { if(values != NULL) { memcpy(this->values, mat.values, m*n*sizeof(Real)); } else { this->m = mat.m; this->n = mat.n; this->values = mat.values; const_cast(mat).wrapped = true; this->wrapped = false; } } return *this; }; /* ---------------------------------------------------------------------- */ inline Matrix & operator=(Real x) { Real * a = this->values; for (UInt i = 0; i < n*m; ++i) *(a++) = x; return *this; } /* ---------------------------------------------------------------------- */ inline Matrix operator* (const Matrix & B) { Matrix C(this->m, B.n); C.mul(*this, B); return C; }; /* ---------------------------------------------------------------------- */ inline Matrix & operator+=(const Matrix & A) { Real * a = this->storage(); Real * b = A.storage(); for (UInt i = 0; i < n*m; ++i) *(a++) += *(b++); return *this; } /* ---------------------------------------------------------------------- */ inline Matrix & operator+=(Real x) { Real * a = this->values; for (UInt i = 0; i < n*m; ++i) *(a++) += x; return *this; } /* ---------------------------------------------------------------------- */ inline Matrix & operator*=(Real x) { Real * a = this->storage(); for (UInt i = 0; i < n*m; ++i) *(a++) *= x; return *this; } /* ---------------------------------------------------------------------- */ inline Matrix & operator/=(Real x) { Real * a = this->values; for (UInt i = 0; i < n*m; ++i) *(a++) /= x; return *this; } /* ---------------------------------------------------------------------- */ template inline void mul(const Matrix & A, const Matrix & B, Real alpha = 1.0) { UInt k = A.n; if(tr_A) k = A.m; #ifndef AKANTU_NDEBUG if (tr_B){ AKANTU_DEBUG_ASSERT(k == B.n, "matrices to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(n == B.m, "matrices to multiply have no fit dimensions"); } else { AKANTU_DEBUG_ASSERT(k == B.m, "matrices to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(n == B.n, "matrices to multiply have no fit dimensions"); } if (tr_A){ AKANTU_DEBUG_ASSERT(m == A.n, "matrices to multiply have no fit dimensions"); } else{ AKANTU_DEBUG_ASSERT(m == A.m, "matrices to multiply have no fit dimensions"); } #endif //AKANTU_NDEBUG Math::matMul(m, n, k, alpha, A.storage(), B.storage(), 0., values); } /* ---------------------------------------------------------------------- */ inline void outerProduct(const types::Vector & A, const types::Vector & B) { AKANTU_DEBUG_ASSERT(A.size() == m && B.size() == n, "A and B are not compatible with the size of the matrix"); for (UInt i = 0; i < m; ++i) { for (UInt j = 0; j < n; ++j) { values[i * m + j] += A[i] * B[j]; } } } /* ---------------------------------------------------------------------- */ inline void eig(types::Vector & eigenvalues, Matrix & eigenvectors) const { AKANTU_DEBUG_ASSERT(n == m, "eig is not a valid operation on a rectangular matrix"); Math::matrixEig(this->n, this->values, eigenvalues.storage(), eigenvectors.storage()); } /* ---------------------------------------------------------------------- */ inline void clear() { std::fill_n(values, m * n, 0); }; /* ---------------------------------------------------------------------- */ inline void copy(const Matrix & src) { memcpy(values, src.storage(), m * n * sizeof(Real)); } /* ---------------------------------------------------------------------- */ inline void eye(Real alpha = 1.) { AKANTU_DEBUG_ASSERT(n == m, "eye is not a valid operation on a rectangular matrix"); clear(); for (UInt i = 0; i < n; ++i) { values[i*n + i] = alpha; } } /* ---------------------------------------------------------------------- */ inline Real trace() const { AKANTU_DEBUG_ASSERT(n == m, "trace is not a valid operation on a rectangular matrix"); Real trace = 0.; for (UInt i = 0; i < n; ++i) trace += values[i*n + i]; return trace; } /* -------------------------------------------------------------------------- */ inline Real norm() const { return Math::norm(this->n*this->m, this->values); } /* ---------------------------------------------------------------------- */ inline Matrix transpose() const { Matrix tmp(m, n); for (UInt i = 0; i < n; ++i) { for (UInt j = 0; j < m; ++j) { tmp(j,i) = operator()(i,j); } } return tmp; } /* ---------------------------------------------------------------------- */ - inline void inv(const Matrix & A) { + inline void inverse(const Matrix & A) { AKANTU_DEBUG_ASSERT(A.n == A.m, "inv is not a valid operation on a rectangular matrix"); AKANTU_DEBUG_ASSERT(n == A.n, "the matrix should have the same size as its inverse"); Math::inv(A.n, A.values, this->values); } /* ---------------------------------------------------------------------- */ /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "types::Matrix" << " [" << n << "," << m <<"] :" << std::endl; for (UInt i = 0; i < m; ++i) { stream << space << AKANTU_INDENT << "| "; for (UInt j = 0; j < n; ++j) { stream << std::setw(10) << values[i*n +j] << " "; } stream << "|" << std::endl; } }; friend class ::akantu::Vector; protected: UInt m; UInt n; Real* values; bool wrapped; }; /* ------------------------------------------------------------------------ */ template template inline void Vector::mul(const Matrix & A, const Vector & x, Real alpha) { UInt n = x.n; Math::matVectMul(this->n, n, alpha, A.storage(), x.storage(), 0., values); } /* -------------------------------------------------------------------------- */ inline std::ostream & operator<<(std::ostream & stream, const Matrix & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Vector & _this) { _this.printself(stream); return stream; } } __END_AKANTU__ #endif /* __AKANTU_AKA_TYPES_HH__ */ diff --git a/src/common/aka_vector.cc b/src/common/aka_vector.cc index c1e341bc1..267f0bd35 100644 --- a/src/common/aka_vector.cc +++ b/src/common/aka_vector.cc @@ -1,139 +1,121 @@ /** * @file aka_vector.cc * @author Nicolas Richart * @date Thu Jun 17 15:14:24 2010 * * @brief class vector * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_vector.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Functions VectorBase */ /* -------------------------------------------------------------------------- */ VectorBase::VectorBase(const ID & id) : id(id), allocated_size(0), size(0), nb_component(1), size_of_type(0) { } /* -------------------------------------------------------------------------- */ VectorBase::~VectorBase() { } /* -------------------------------------------------------------------------- */ void VectorBase::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "VectorBase [" << std::endl; stream << space << " + size : " << size << std::endl; stream << space << " + nb component : " << nb_component << std::endl; stream << space << " + allocated size : " << allocated_size << std::endl; Real mem_size = (allocated_size * nb_component * size_of_type) / 1024.; stream << space << " + size of type : " << size_of_type << "B" << std::endl; stream << space << " + memory allocated : " << mem_size << "kB" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template <> Int Vector::find(const Real & elem) const { AKANTU_DEBUG_IN(); UInt i = 0; Real epsilon = std::numeric_limits::epsilon(); for (; (i < size) && (fabs(values[i] - elem) <= epsilon); ++i); AKANTU_DEBUG_OUT(); return (i == size) ? -1 : (Int) i; } /* -------------------------------------------------------------------------- */ template <> Vector & Vector::operator*=(__attribute__((unused)) const ElementType & alpha) { AKANTU_DEBUG_TO_IMPLEMENT(); return *this; } template <> Vector & Vector::operator-=(__attribute__((unused)) const Vector & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); return *this; } template <> Vector & Vector::operator+=(__attribute__((unused)) const Vector & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); return *this; } template <> Vector & Vector::operator*=(__attribute__((unused)) const char & alpha) { AKANTU_DEBUG_TO_IMPLEMENT(); return *this; } template <> Vector & Vector::operator-=(__attribute__((unused)) const Vector & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); return *this; } template <> Vector & Vector::operator+=(__attribute__((unused)) const Vector & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); return *this; } - -/* -------------------------------------------------------------------------- */ -#define AKANTU_DESTRUCTOR_SPEC(type) \ - template <> Vector::~Vector () { \ - AKANTU_DEBUG_IN(); \ - AKANTU_DEBUG(dblAccessory, "Freeing " \ - << allocated_size*nb_component*sizeof(type) / 1024. \ - << "kB (" << id <<")"); \ - if(values) free(values); \ - size = allocated_size = 0; \ - AKANTU_DEBUG_OUT(); \ - } -AKANTU_DESTRUCTOR_SPEC(Real) -AKANTU_DESTRUCTOR_SPEC(UInt) -AKANTU_DESTRUCTOR_SPEC(Int) -AKANTU_DESTRUCTOR_SPEC(bool) - - /* -------------------------------------------------------------------------- */ // template class Vector; // template class Vector; // template class Vector; // template class Vector; // template class Vector; // template class Vector; // template class Vector; //template class Vector; __END_AKANTU__ diff --git a/src/common/aka_vector.hh b/src/common/aka_vector.hh index dd92371db..2d60d0cb5 100644 --- a/src/common/aka_vector.hh +++ b/src/common/aka_vector.hh @@ -1,356 +1,355 @@ /** * @file aka_vector.hh * @author Nicolas Richart * @date Thu Jun 17 10:04:55 2010 * * @brief class of vectors * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_VECTOR_HH__ #define __AKANTU_VECTOR_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class Matrix; /// class that afford to store vectors in static memory class VectorBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: VectorBase(const ID & id = ""); virtual ~VectorBase(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get the amount of space allocated in bytes inline UInt getMemorySize() const; /// set the size to zero without freeing the allocated space inline void empty(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(AllocatedSize, allocated_size, UInt); AKANTU_GET_MACRO(Size, size, UInt); AKANTU_GET_MACRO(NbComponent, nb_component, UInt); AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO(Tag, tag, const std::string &); AKANTU_SET_MACRO(Tag, tag, const std::string &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the vector ID id; /// the size allocated UInt allocated_size; /// the size used UInt size; /// number of components UInt nb_component; /// size of the stored type UInt size_of_type; /// User defined tag std::string tag; }; namespace types { class Matrix; template class Vector; } /* -------------------------------------------------------------------------- */ template class Vector : public VectorBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef T value_type; typedef value_type & reference; typedef value_type * pointer_type; typedef const value_type & const_reference; /// Allocation of a new vector - Vector(UInt size = 0, UInt nb_component = 1, - const ID & id = ""); + inline Vector(UInt size = 0, UInt nb_component = 1, + const ID & id = ""); /// Allocation of a new vector with a default value Vector(UInt size, UInt nb_component, const value_type def_values[], const ID & id = ""); /// Allocation of a new vector with a default value Vector(UInt size, UInt nb_component, const_reference value, const ID & id = ""); /// Copy constructor (deep copy if deep=true) Vector(const Vector& vect, bool deep = true, const ID & id = ""); /// Copy constructor (deep copy) Vector(const std::vector& vect); - virtual ~Vector(); + virtual inline ~Vector(); /* ------------------------------------------------------------------------ */ /* Iterator */ /* ------------------------------------------------------------------------ */ /// \todo protected: does not compile with intel check why public: template class iterator_internal; public: /* ------------------------------------------------------------------------ */ // template class iterator : public iterator_internal {}; // template class const_iterator : public iterator_internal {}; /* ------------------------------------------------------------------------ */ template class iterator : public iterator_internal { public: typedef iterator_internal parent; typedef typename parent::pointer pointer; public: iterator() : parent() {}; iterator(pointer_type data, UInt offset) : parent(data, offset) {}; iterator(pointer warped) : parent(warped) {}; iterator(const iterator & it) : parent(it) {}; }; /* ------------------------------------------------------------------------ */ template class const_iterator : public iterator_internal { public: typedef iterator_internal parent; typedef typename parent::pointer pointer; public: const_iterator() : parent() {}; const_iterator(pointer_type data, UInt offset) : parent(data, offset) {}; const_iterator(pointer warped) : parent(warped) {}; const_iterator(const const_iterator & it) : parent(it) {}; }; /// specialization of the previous iterators for the scalar types /* ------------------------------------------------------------------------ */ template class iterator : public iterator_internal { public: typedef iterator_internal parent; typedef typename parent::pointer pointer; public: iterator() : parent() {}; iterator(pointer_type data, UInt offset) : parent(data, offset) {}; iterator(pointer warped) : parent(warped) {}; iterator(const iterator & it) : parent(it) {}; }; /* -------------------------------------------------------------------------- */ template class const_iterator : public iterator_internal { public: typedef iterator_internal parent; typedef typename parent::pointer pointer; public: const_iterator() : parent() {}; const_iterator(pointer_type data, UInt offset) : parent(data, offset) {}; const_iterator(pointer warped) : parent(warped) {}; const_iterator(const const_iterator & it) : parent(it) {}; }; /* ------------------------------------------------------------------------ */ // template inline iterator begin(); // template inline iterator end(); // template inline const_iterator begin() const; // template inline const_iterator end() const; inline iterator begin(); inline iterator end(); inline const_iterator begin() const; inline const_iterator end() const; inline iterator< types::Vector > begin(UInt n); inline iterator< types::Vector > end(UInt n); inline const_iterator< types::Vector > begin(UInt n) const; inline const_iterator< types::Vector > end(UInt n) const; inline iterator< types::Matrix > begin(UInt m, UInt n); inline iterator< types::Matrix > end(UInt m, UInt n); inline const_iterator< types::Matrix > begin(UInt m, UInt n) const; inline const_iterator< types::Matrix > end(UInt m, UInt n) const; inline iterator< types::Matrix > begin_reinterpret(UInt m, UInt n, UInt size, UInt nb_component); inline iterator< types::Matrix > end_reinterpret(UInt m, UInt n, UInt size, UInt nb_component); inline const_iterator< types::Matrix > begin_reinterpret(UInt m, UInt n, UInt size, UInt nb_component) const; inline const_iterator< types::Matrix > end_reinterpret(UInt m, UInt n, UInt size, UInt nb_component) const; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get jth componemt of the ith tuple in read-only inline const_reference get(UInt i, UInt j = 0) const; /// get jth component of the ith tuple inline reference at(UInt i, UInt j = 0); /// add an element at the end of the vector with the value value for all /// component inline void push_back(const_reference value); /// add an element at the end of the vector inline void push_back(const value_type new_elem[]); /** * remove an element and move the last one in the hole * /!\ change the order in the vector */ inline void erase(UInt i); template inline void erase(const iterator & it); /// change the size of the vector and allocate more memory if needed void resize(UInt size); /// change the number of components by interlacing data void extendComponentsInterlaced(UInt multiplicator, UInt stride); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; // Vector& operator=(const Vector& vect); /// search elem in the vector, return the position of the first occurrence or /// -1 if not found Int find(const_reference elem) const; Int find(T elem[]) const; /// set a vvector to 0 inline void clear() { std::fill_n(values, size*nb_component, T()); }; /// copy the content of an other vector void copy(const Vector & vect); /// give the address of the memory allocated for this vector T * storage() const { return values; }; protected: /// perform the allocation for the constructors void allocate(UInt size, UInt nb_component = 1); /// resize without initializing the memory void resizeUnitialized(UInt new_size); /* ------------------------------------------------------------------------ */ /* Operators */ /* ------------------------------------------------------------------------ */ public: Vector & operator-=(const Vector & vect); Vector & operator+=(const Vector & vect); Vector & operator*=(const T & alpha); inline reference operator()(UInt i, UInt j = 0); inline const_reference operator()(UInt i, UInt j = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: UInt getSize() const{ return this->size; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: /// array of values T * values; // /!\ very dangerous }; __END_AKANTU__ #include "aka_types.hh" - __BEGIN_AKANTU__ #include "aka_vector_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Inline Functions Vector */ /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Vector & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ /* Inline Functions VectorBase */ /* -------------------------------------------------------------------------- */ inline std::ostream & operator<<(std::ostream & stream, const VectorBase & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_VECTOR_HH__ */ diff --git a/src/common/aka_vector_tmpl.hh b/src/common/aka_vector_tmpl.hh index f3e89b5cd..09ec963be 100644 --- a/src/common/aka_vector_tmpl.hh +++ b/src/common/aka_vector_tmpl.hh @@ -1,1051 +1,1085 @@ /** * @file aka_vector_inline_impl.cc * @author Nicolas Richart * @date Thu Jul 15 00:09:33 2010 * * @brief Inline functions of the classes Vector and VectorBase * * @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 . * */ /* -------------------------------------------------------------------------- */ /* Inline Functions Vector */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template inline T & Vector::operator()(UInt i, UInt j) { AKANTU_DEBUG_ASSERT(size > 0, "The vector \"" << id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range in vector \"" << id << "\""); return values[i*nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline const T & Vector::operator()(UInt i, UInt j) const { AKANTU_DEBUG_ASSERT(size > 0, "The vector \"" << id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range in vector \"" << id << "\""); return values[i*nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline T & Vector::at(UInt i, UInt j) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(size > 0, "The vector is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range"); AKANTU_DEBUG_OUT(); return values[i*nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline const T & Vector::get(UInt i, UInt j) const{ AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(size > 0, "The vector is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range"); AKANTU_DEBUG_OUT(); return values[i*nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline void Vector::push_back(const T & value) { // AKANTU_DEBUG_IN(); UInt pos = size; resizeUnitialized(size+1); /// @todo see if with std::uninitialized_fill it allow to build vector of objects std::uninitialized_fill_n(values + pos * nb_component, nb_component, value); // for (UInt i = 0; i < nb_component; ++i) { // values[pos*nb_component + i] = value; // } // AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void Vector::push_back(const T new_elem[]) { // AKANTU_DEBUG_IN(); UInt pos = size; resizeUnitialized(size+1); T * tmp = values + nb_component * pos; std::uninitialized_copy(new_elem, new_elem + nb_component, tmp); // for (UInt i = 0; i < nb_component; ++i) { // values[pos*nb_component + i] = new_elem[i]; // } // AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void Vector::erase(UInt i){ AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT((size > 0), "The vector is empty"); AKANTU_DEBUG_ASSERT((i < size), "The element at position [" << i << "] is out of range"); if(i != (size - 1)) { for (UInt j = 0; j < nb_component; ++j) { values[i*nb_component + j] = values[(size-1)*nb_component + j]; } } resize(size - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector & Vector::operator-=(const Vector & vect) { AKANTU_DEBUG_ASSERT((size == vect.size) && (nb_component == vect.nb_component), "The too vector don't have the same sizes"); T * a = values; T * b = vect.values; for (UInt i = 0; i < size*nb_component; ++i) { *a -= *b; ++a;++b; } return *this; } /* -------------------------------------------------------------------------- */ template Vector & Vector::operator+=(const Vector & vect) { AKANTU_DEBUG_ASSERT((size == vect.size) && (nb_component == vect.nb_component), "The too vector don't have the same sizes"); T * a = values; T * b = vect.values; for (UInt i = 0; i < size*nb_component; ++i) { *a++ += *b++; } return *this; } /* -------------------------------------------------------------------------- */ template Vector & Vector::operator*=(const T & alpha) { T * a = values; for (UInt i = 0; i < size*nb_component; ++i) { *a++ *= alpha; } return *this; } /* -------------------------------------------------------------------------- */ /* Functions Vector */ /* -------------------------------------------------------------------------- */ -template Vector::Vector (UInt size, - UInt nb_component, - const ID & id) : +template inline Vector::Vector (UInt size, + UInt nb_component, + const ID & id) : VectorBase(id), values(NULL) { AKANTU_DEBUG_IN(); allocate(size, nb_component); T val = T(); std::uninitialized_fill(values, values + size*nb_component, val); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +#define AKANTU_CONSTRUCTOR_SPEC(type) \ + template<> inline Vector::Vector (UInt size, \ + UInt nb_component, \ + const ID & id) : \ + VectorBase(id), values(NULL) { \ + AKANTU_DEBUG_IN(); \ + allocate(size, nb_component); \ + AKANTU_DEBUG_OUT(); \ + } + +AKANTU_CONSTRUCTOR_SPEC(Real) +AKANTU_CONSTRUCTOR_SPEC(UInt) +AKANTU_CONSTRUCTOR_SPEC(Int) +AKANTU_CONSTRUCTOR_SPEC(bool) +#undef AKANTU_CONSTRUCTOR_SPEC + + /* -------------------------------------------------------------------------- */ template Vector::Vector (UInt size, UInt nb_component, const T def_values[], const ID & id) : VectorBase(id), values(NULL) { AKANTU_DEBUG_IN(); allocate(size, nb_component); T * tmp = values; for (UInt i = 0; i < size; ++i) { tmp = values + nb_component * i; std::uninitialized_copy(def_values, def_values + nb_component, tmp); // for (UInt j = 0; j < nb_component; ++j) { // values[i*nb_component + j] = def_values[j]; // } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector::Vector (UInt size, UInt nb_component, const T & value, const ID & id) : VectorBase(id), values(NULL) { AKANTU_DEBUG_IN(); allocate(size, nb_component); std::uninitialized_fill_n(values, size*nb_component, value); // for (UInt i = 0; i < nb_component*size; ++i) { // values[i] = value; // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector::Vector(const Vector& vect, bool deep, const ID & id) { AKANTU_DEBUG_IN(); this->id = (id == "") ? vect.id : id; if (deep) { allocate(vect.size, vect.nb_component); T * tmp = values; std::uninitialized_copy(vect.values, vect.values + size * nb_component, tmp); // for (UInt i = 0; i < size; ++i) { // for (UInt j = 0; j < nb_component; ++j) { // values[i*nb_component + j] = vect.values[i*nb_component + j]; // } // } } else { this->values = vect.values; this->size = vect.size; this->nb_component = vect.nb_component; this->allocated_size = vect.allocated_size; this->size_of_type = vect.size_of_type; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector::Vector(const std::vector& vect) { AKANTU_DEBUG_IN(); this->id = ""; allocate(vect.size(), 1); T * tmp = values; std::uninitialized_copy(&(vect[0]), &(vect[size-1]), tmp); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template Vector::~Vector () { +template inline Vector::~Vector () { AKANTU_DEBUG_IN(); AKANTU_DEBUG(dblAccessory, "Freeing " << allocated_size*nb_component*sizeof(T) / 1024. << "kB (" << id <<")"); if(values){ for (UInt i = 0; i < size * nb_component; ++i) { T * obj = values+i; obj->~T(); } free(values); } size = allocated_size = 0; AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +#define AKANTU_DESTRUCTOR_SPEC(type) \ + template<> inline Vector::~Vector () { \ + AKANTU_DEBUG_IN(); \ + AKANTU_DEBUG(dblAccessory, "Freeing " \ + << allocated_size*nb_component*sizeof(type) / 1024. \ + << "kB (" << id <<")"); \ + if(values) free(values); \ + size = allocated_size = 0; \ + AKANTU_DEBUG_OUT(); \ + } +AKANTU_DESTRUCTOR_SPEC(Real) +AKANTU_DESTRUCTOR_SPEC(UInt) +AKANTU_DESTRUCTOR_SPEC(Int) +AKANTU_DESTRUCTOR_SPEC(bool) +#undef AKANTU_DESTRUCTOR_SPEC /* -------------------------------------------------------------------------- */ template void Vector::allocate(UInt size, UInt nb_component) { AKANTU_DEBUG_IN(); if (size == 0){ values = NULL; } else { values = static_cast(malloc(nb_component * size * sizeof(T))); AKANTU_DEBUG_ASSERT(values != NULL, "Cannot allocate " << nb_component * size * sizeof(T) / 1024. << "kB (" << id <<")"); } if (values == NULL) { this->size = this->allocated_size = 0; } else { AKANTU_DEBUG(dblAccessory, "Allocated " << size * nb_component * sizeof(T) / 1024. << "kB (" << id <<")"); this->size = this->allocated_size = size; } this->size_of_type = sizeof(T); this->nb_component = nb_component; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Vector::resize(UInt new_size) { UInt old_size = size; T * old_values = values; if(new_size < size) { for (UInt i = new_size * nb_component; i < size * nb_component; ++i) { T * obj = old_values+i; obj->~T(); } } resizeUnitialized(new_size); T val = T(); if(size > old_size) std::uninitialized_fill(values + old_size*nb_component, values + size*nb_component, val); } /* -------------------------------------------------------------------------- */ template void Vector::resizeUnitialized(UInt new_size) { // AKANTU_DEBUG_IN(); // free some memory if(new_size <= allocated_size) { if(allocated_size - new_size > AKANTU_MIN_ALLOCATION) { AKANTU_DEBUG(dblAccessory, "Freeing " << (allocated_size - size)*nb_component*sizeof(T) / 1024. << "kB (" << id <<")"); // Normally there are no allocation problem when reducing an array T * tmp_ptr = static_cast(realloc(values, new_size * nb_component * sizeof(T))); if(new_size != 0 && tmp_ptr == NULL) { AKANTU_DEBUG_ERROR("Cannot free data (" << id << ")" << " [current allocated size : " << allocated_size << " | " << "requested size : " << new_size << "]"); } values = tmp_ptr; allocated_size = new_size; } size = new_size; // AKANTU_DEBUG_OUT(); return; } // allocate more memory UInt size_to_alloc = (new_size - allocated_size < AKANTU_MIN_ALLOCATION) ? allocated_size + AKANTU_MIN_ALLOCATION : new_size; T *tmp_ptr = static_cast(realloc(values, size_to_alloc * nb_component * sizeof(T))); AKANTU_DEBUG_ASSERT(tmp_ptr != NULL, "Cannot allocate " << size_to_alloc * nb_component * sizeof(T) / 1024. << "kB"); if (tmp_ptr == NULL) { AKANTU_DEBUG_ERROR("Cannot allocate more data (" << id << ")" << " [current allocated size : " << allocated_size << " | " << "requested size : " << new_size << "]"); } AKANTU_DEBUG(dblAccessory, "Allocating " << (size_to_alloc - allocated_size)*nb_component*sizeof(T) / 1024. << "kB"); allocated_size = size_to_alloc; size = new_size; values = tmp_ptr; // AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Vector::extendComponentsInterlaced(UInt multiplicator, UInt block_size) { AKANTU_DEBUG_IN(); if (multiplicator == 1) return; AKANTU_DEBUG_ASSERT(multiplicator > 1, "invalid multiplicator"); AKANTU_DEBUG_ASSERT(nb_component%block_size == 0, "stride must divide actual number of components"); values = static_cast(realloc(values, nb_component*multiplicator*size* sizeof(T))); UInt new_component = nb_component/block_size * multiplicator; for (UInt i = 0,k=size-1; i < size; ++i,--k) { for (UInt j = 0; j < new_component; ++j) { UInt m = new_component - j -1; UInt n = m/multiplicator; for (UInt l = 0,p=block_size-1; l < block_size; ++l,--p) { values[k*nb_component*multiplicator+m*block_size+p] = values[k*nb_component+n*block_size+p]; } } } nb_component = nb_component * multiplicator; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Int Vector::find(const T & elem) const { AKANTU_DEBUG_IN(); UInt i = 0; for (; (i < size) && (values[i] != elem); ++i); AKANTU_DEBUG_OUT(); return (i == size) ? -1 : (Int) i; } /* -------------------------------------------------------------------------- */ template Int Vector::find(T elem[]) const { AKANTU_DEBUG_IN(); T * it = values; UInt i = 0; UInt c = 0; for (;i < size && (c != nb_component); ++i) { c = 0; T * cit = it; T * celem = elem; for(; (c < nb_component) && (*cit == *celem); ++c, ++cit, ++celem); it += nb_component; } AKANTU_DEBUG_OUT(); return (i == size) ? -1 : (Int) i; } /* -------------------------------------------------------------------------- */ template void Vector::copy(const Vector& vect) { AKANTU_DEBUG_IN(); if(AKANTU_DEBUG_TEST(dblWarning)) if(vect.nb_component != nb_component) { AKANTU_DEBUG(dblWarning, "The two vectors do not have the same number of components"); } // this->id = vect.id; resize((vect.size * vect.nb_component) / nb_component); T * tmp = values; std::uninitialized_copy(vect.values, vect.values + size * nb_component, tmp); // memcpy(this->values, vect.values, vect.size * vect.nb_component * sizeof(T)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Vector::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); Real real_size = allocated_size * nb_component * size_of_type / 1024.0; std::streamsize prec = stream.precision(); std::ios_base::fmtflags ff = stream.flags(); stream.setf (std::ios_base::showbase); stream.precision(2); stream << space << "Vector<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; stream << space << " + id : " << this->id << std::endl; stream << space << " + size : " << this->size << std::endl; stream << space << " + nb_component : " << this->nb_component << std::endl; stream << space << " + allocated size : " << this->allocated_size << std::endl; stream << space << " + memory size : " << real_size << "kB" << std::endl; if(!AKANTU_DEBUG_LEVEL_IS_TEST()) stream << space << " + address : " << std::hex << this->values << std::dec << std::endl; stream.precision(prec); stream.flags(ff); if(AKANTU_DEBUG_TEST(dblDump)) { stream << space << " + values : {"; for (UInt i = 0; i < this->size; ++i) { stream << "{"; for (UInt j = 0; j < this->nb_component; ++j) { stream << this->values[i*nb_component + j]; if(j != this->nb_component - 1) stream << ", "; } stream << "}"; if(i != this->size - 1) stream << ", "; } stream << "}" << std::endl; } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /* Inline Functions VectorBase */ /* -------------------------------------------------------------------------- */ inline UInt VectorBase::getMemorySize() const { return allocated_size * nb_component * size_of_type; } inline void VectorBase::empty() { size = 0; } /** \todo change definition of iterators template template template template::isScalar > class Vector::iterator_internal : public std::iterator { */ /* -------------------------------------------------------------------------- */ /* Iterators */ /* -------------------------------------------------------------------------- */ template template //class Vector::iterator_internal : public virtual std::iterator { class Vector::iterator_internal { public: typedef R value_type; typedef R* pointer; typedef R& reference; typedef IR internal_value_type; typedef IR* internal_pointer; protected: iterator_internal(UInt offset, pointer_type data, pointer ret) : offset(offset), initial(data), ret(ret) { } virtual ~iterator_internal() { delete ret; }; public: iterator_internal() : offset(0), initial(NULL), ret(NULL) {}; iterator_internal(pointer_type data, UInt offset) : offset(offset), initial(data), ret(new value_type(data)) { AKANTU_DEBUG_ASSERT(offset == ret->size(), "The iterator_internal is not compatible with the type " << typeid(value_type).name()); } iterator_internal(pointer warped) : offset(warped->size()), initial(warped->storage()), ret(const_cast(warped)) { } iterator_internal(const iterator_internal & it) { if(this != &it) { this->offset = it.offset; this->initial = it.initial; this->ret = new internal_value_type(*it.ret); } } inline iterator_internal & operator=(const iterator_internal & it) { if(this != &it) { this->offset = it.offset; this->initial = it.initial; if(this->ret) this->ret->shallowCopy(*it.ret); else this->ret = new internal_value_type(*it.ret); } return *this; } inline reference operator*() { return *ret; }; inline pointer operator->() { return ret; }; inline iterator_internal & operator++() { ret->values += offset; return *this; }; inline iterator_internal & operator+=(const UInt n) { ret->values += offset * n; return *this; } inline reference operator[](const UInt n) { ret->values = initial + n*offset; return *ret; } inline bool operator==(const iterator_internal & other) { return (*this).ret->values == other.ret->values; } inline bool operator!=(const iterator_internal & other) { return (*this).ret->values != other.ret->values; } inline pointer_type getCurrentStorage() const { return ret->storage(); } protected: UInt offset; pointer_type initial; internal_pointer ret; }; // /* -------------------------------------------------------------------------- */ // template // template // inline Vector::iterator Vector::begin() { // return iterator(values, nb_component); // } // /* -------------------------------------------------------------------------- */ // template // template // inline Vector::iterator Vector::end() { // return iterator(values + nb_component * size, nb_component); // } /* -------------------------------------------------------------------------- */ template inline Vector::iterator< types::Vector > Vector::begin(UInt n) { AKANTU_DEBUG_ASSERT(nb_component == n, "The iterator is not compatible with the type Vector(" << n<< ")"); return iterator< types::Vector >(new types::Vector(values, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::iterator< types::Vector > Vector::end(UInt n) { AKANTU_DEBUG_ASSERT(nb_component == n, "The iterator is not compatible with the type Vector(" << n<< ")"); return iterator< types::Vector >(new types::Vector(values + nb_component * size, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator< types::Vector > Vector::begin(UInt n) const { AKANTU_DEBUG_ASSERT(nb_component == n, "The iterator is not compatible with the type Vector(" << n<< ")"); return const_iterator< types::Vector >(new types::Vector(values, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator< types::Vector > Vector::end(UInt n) const { AKANTU_DEBUG_ASSERT(nb_component == n, "The iterator is not compatible with the type Vector(" << n<< ")"); return const_iterator< types::Vector >(new types::Vector(values + nb_component * size, n)); } /* -------------------------------------------------------------------------- */ template<> inline Vector::iterator Vector::begin(UInt m, UInt n) { AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return iterator< types::Matrix >(new types::Matrix(values, m, n)); } /* -------------------------------------------------------------------------- */ template<> inline Vector::iterator Vector::end(UInt m, UInt n) { AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return iterator(new types::Matrix(values + nb_component * size, m, n)); } /* -------------------------------------------------------------------------- */ template<> inline Vector::const_iterator Vector::begin(UInt m, UInt n) const { AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return const_iterator< types::Matrix >(new types::Matrix(values, m, n)); } /* -------------------------------------------------------------------------- */ template<> inline Vector::const_iterator Vector::end(UInt m, UInt n) const { AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return const_iterator(new types::Matrix(values + nb_component * size, m, n)); } /* -------------------------------------------------------------------------- */ template<> inline Vector::iterator< types::Matrix > Vector::begin_reinterpret(UInt m, UInt n, __attribute__((unused)) UInt size, __attribute__((unused)) UInt nb_component) { AKANTU_DEBUG_ASSERT(nb_component * size == this->nb_component * this->size, "The new values for size (" << size << ") and nb_component (" << nb_component << ") are not compatible with the one of this vector"); AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return iterator< types::Matrix >(new types::Matrix(values + nb_component * 0, m, n)); } /* -------------------------------------------------------------------------- */ template<> inline Vector::iterator< types::Matrix > Vector::end_reinterpret(UInt m, UInt n, UInt size, UInt nb_component) { AKANTU_DEBUG_ASSERT(nb_component * size == this->nb_component * this->size, "The new values for size (" << size << ") and nb_component (" << nb_component << ") are not compatible with the one of this vector"); AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return iterator(new types::Matrix(values + nb_component * size, m, n)); } /* -------------------------------------------------------------------------- */ template<> inline Vector::const_iterator< types::Matrix > Vector::begin_reinterpret(UInt m, UInt n, __attribute__((unused)) UInt size, __attribute__((unused)) UInt nb_component) const { AKANTU_DEBUG_ASSERT(nb_component * size == this->nb_component * this->size, "The new values for size (" << size << ") and nb_component (" << nb_component << ") are not compatible with the one of this vector"); AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return const_iterator< types::Matrix >(new types::Matrix(values + nb_component * 0, m, n)); } /* -------------------------------------------------------------------------- */ template<> inline Vector::const_iterator< types::Matrix > Vector::end_reinterpret(UInt m, UInt n, UInt size, UInt nb_component) const { AKANTU_DEBUG_ASSERT(nb_component * size == this->nb_component * this->size, "The new values for size (" << size << ") and nb_component (" << nb_component << ") are not compatible with the one of this vector"); AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return const_iterator(new types::Matrix(values + nb_component * size, m, n)); } // template // template // inline typename Vector::template iterator_internal & Vector::iterator_internal::operator++() { // ret += offset; // return *this; // } // inline iterator_internal & operator+=(const UInt n) { // ret->values += offset * n; // return *this; // } /* -------------------------------------------------------------------------- */ /** * Specialization for scalar types */ /* If gcc < 4.4 some problem with detection of specialization of a template by * an other template so I have copied the code, that is really ugly just to help * the compiler. This part of code should be removed if version of gcc before * * 4.4 are note supported anymore. Problems for know appears on Mac OS with gcc * version is 4.2.1 */ #if (__GNUC__ < 4) || \ ((__GNUC__ == 4) && (__GNUC_MINOR__ < 4)) #define specialize_internal_iterator_for_scalar(T) \ template <> \ template <> \ class Vector::iterator_internal { \ public: \ typedef T value_type; \ typedef T* pointer; \ typedef T& reference; \ typedef T internal_value_type; \ typedef T* internal_pointer; \ protected: \ iterator_internal(__attribute__((unused)) UInt offset, \ pointer data, pointer ret) : \ initial(data), \ ret(ret) { \ AKANTU_DEBUG_ASSERT(offset == 1, \ "The iterator is not compatible with the type " \ << typeid(value_type).name()); \ } \ \ public: \ iterator_internal() : initial(NULL), ret(NULL) {}; \ \ iterator_internal(pointer data, \ __attribute__((unused)) UInt offset) : \ initial(data), \ ret(data) { \ AKANTU_DEBUG_ASSERT(offset == 1, \ "The iterator is not compatible with the type " \ << typeid(value_type).name()); \ }; \ \ iterator_internal(const iterator_internal & it) { \ if(this != &it) { \ this->initial = it.initial; \ this->ret = new internal_value_type(*it.ret); \ } \ } \ \ virtual ~iterator_internal() { }; \ \ inline iterator_internal & operator=(const iterator_internal & it) { \ if(this != &it) { \ this->initial = it.initial; \ this->ret = it.ret; \ } \ return *this; \ } \ \ inline reference operator*() { return *ret; }; \ inline pointer operator->() { return ret; }; \ inline iterator_internal & operator++() { ++ret; return *this; }; \ \ inline iterator_internal & operator+=(const UInt n) { \ ret += n; \ return *this; \ } \ \ inline reference operator[](const UInt n) { \ ret = initial + n; \ return *ret; \ } \ \ inline bool operator==(const iterator_internal & other) { \ return (*this).ret == other.ret; \ } \ inline bool operator!=(const iterator_internal & other) { \ return (*this).ret != other.ret; \ } \ \ inline iterator_internal & operator-(size_t n) { \ ret -= n; \ return *this; \ } \ \ inline size_t operator-(const iterator_internal & b) { \ return ret - b.getCurrentStorage(); \ } \ \ inline pointer getCurrentStorage() const { \ return ret; \ } \ \ protected: \ pointer initial; \ pointer ret; \ } specialize_internal_iterator_for_scalar(Real); specialize_internal_iterator_for_scalar(UInt); #else template template class Vector::iterator_internal : public std::iterator{ public: typedef T value_type; typedef T* pointer; typedef T& reference; typedef T internal_value_type; typedef T* internal_pointer; protected: iterator_internal(UInt offset, pointer data, pointer ret) : initial(data), ret(ret) { } public: iterator_internal() : initial(NULL), ret(NULL) {}; iterator_internal(pointer data, UInt offset) : initial(data), ret(data) { AKANTU_DEBUG_ASSERT(offset == 1, "The iterator_internal is not compatible with the type " << typeid(value_type).name()); }; iterator_internal(const iterator_internal & it) { if(this != &it) { this->initial = it.initial; this->ret = new internal_value_type(*it.ret); } } virtual ~iterator_internal() { }; inline iterator_internal & operator=(const iterator_internal & it) { if(this != &it) { this->initial = it.initial; this->ret = it.ret; } return *this; } inline reference operator*() { return *ret; }; inline pointer operator->() { return ret; }; inline iterator_internal & operator++() { ++ret; return *this; }; inline iterator_internal & operator--() { --ret; return *this; }; inline iterator_internal & operator+=(const UInt n) { ret += n; return *this; } inline reference operator[](const UInt n) { ret = initial + n; return *ret; } inline bool operator==(const iterator_internal & other) { return (*this).ret == other.ret; } inline bool operator!=(const iterator_internal & other) { return (*this).ret != other.ret; } inline bool operator<(const iterator_internal & other) { return ret < other.ret; } inline bool operator<=(const iterator_internal & other) { return ret <= other.ret; } inline bool operator>(const iterator_internal & other) { return ret > other.ret; } inline bool operator>=(const iterator_internal & other) { return ret >= other.ret; } inline iterator_internal & operator-(size_t n) { ret -= n; return *this; } inline size_t operator-(const iterator_internal & b) { return ret - b.getCurrentStorage(); } inline iterator_internal & operator+(size_t n) { ret += n; return *this; } inline pointer getCurrentStorage() const { return ret; } protected: pointer initial; pointer ret; }; #endif /* -------------------------------------------------------------------------- */ template template inline void Vector::erase(const iterator & it) { T * curr = it.getCurrentStorage(); UInt pos = (curr - values) / nb_component; erase(pos); } // inline reference operator[](const UInt n) { // ret->values = initial + n*offset; // return *ret; // } /* -------------------------------------------------------------------------- */ template inline Vector::iterator Vector::begin() { return iterator(values, 1); } /* -------------------------------------------------------------------------- */ template inline Vector::iterator Vector::end() { return iterator(values + nb_component * size, 1); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator Vector::begin() const { return const_iterator(values, 1); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator Vector::end() const { return const_iterator(values + nb_component * size, 1); } diff --git a/src/fem/mesh.hh b/src/fem/mesh.hh index d2677ddd2..bf052f4c4 100644 --- a/src/fem/mesh.hh +++ b/src/fem/mesh.hh @@ -1,636 +1,637 @@ /** * @file mesh.hh * @author Nicolas Richart * @date Wed Jun 16 11:53:53 2010 * * @brief the class representing the meshes * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_vector.hh" #include "element_class.hh" #include "by_element_type.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Element */ /* -------------------------------------------------------------------------- */ class Element { public: Element(ElementType type = _not_defined, UInt element = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) : type(type), element(element), ghost_type(ghost_type), kind(kind) {}; Element(const Element & element) { this->type = element.type; this->element = element.element; this->ghost_type = element.ghost_type; this->kind = element.kind; } inline bool operator==(const Element & elem) const { return ((element == elem.element) && (type == elem.type) && (ghost_type == elem.ghost_type) && (kind == elem.kind)); } inline bool operator!=(const Element & elem) const { return ((element != elem.element) || (type != elem.type) || (ghost_type != elem.ghost_type) || (kind != elem.kind)); } virtual ~Element() {}; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; public: ElementType type; UInt element; GhostType ghost_type; ElementKind kind; }; struct CompElementLess { bool operator() (const Element& lhs, const Element& rhs) const { bool res = ((lhs.kind < rhs.kind) || ((lhs.kind == rhs.kind) && ((lhs.ghost_type < rhs.ghost_type) || ((lhs.ghost_type == rhs.ghost_type) && ((lhs.type < rhs.type) || ((lhs.type == rhs.type) && (lhs.element < rhs.element))))))); return res; } }; extern const Element ElementNull; /* -------------------------------------------------------------------------- */ /* ByElementType */ /* -------------------------------------------------------------------------- */ template class ByElementType { protected: typedef std::map DataMap; public: ByElementType(const ID & id = "by_element_type", const ID & parent_id = ""); ~ByElementType(); inline static std::string printType(const ElementType & type, const GhostType & ghost_type); inline bool exists(ElementType type, GhostType ghost_type = _not_ghost) const; inline const Stored & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; inline Stored & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost); inline Stored & operator()(const Stored & insert, const ElementType & type, const GhostType & ghost_type = _not_ghost); void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ class type_iterator : private std::iterator { public: typedef const ElementType value_type; typedef const ElementType* pointer; typedef const ElementType& reference; protected: typedef typename ByElementType::DataMap::const_iterator DataMapIterator; public: type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek); type_iterator(const type_iterator & it); inline reference operator*(); inline type_iterator & operator++(); type_iterator operator++(int); inline bool operator==(const type_iterator & other); inline bool operator!=(const type_iterator & other); private: DataMapIterator list_begin; DataMapIterator list_end; UInt dim; ElementKind kind; }; inline type_iterator firstType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; inline type_iterator lastType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; inline void setID(const ID & id) { this->id = id; } protected: inline DataMap & getData(GhostType ghost_type); inline const DataMap & getData(GhostType ghost_type) const; /* -------------------------------------------------------------------------- */ protected: ID id; DataMap data; DataMap ghost_data; }; /* -------------------------------------------------------------------------- */ /* Some typedefs */ /* -------------------------------------------------------------------------- */ template class ByElementTypeVector : public ByElementType *>, protected Memory { protected: typedef typename ByElementType *>::DataMap DataMap; public: ByElementTypeVector() {}; // ByElementTypeVector(const ID & id = "by_element_type_vector", // const MemoryID & memory_id = 0) : // ByElementType *>(id, memory_id) {}; ByElementTypeVector(const ID & id, const ID & parent_id, const MemoryID & memory_id = 0) : ByElementType *>(id, parent_id), Memory(memory_id) {}; inline Vector & alloc(UInt size, UInt nb_component, const ElementType & type, const GhostType & ghost_type); inline void alloc(UInt size, UInt nb_component, const ElementType & type); inline const Vector & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; inline Vector & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost); inline void setVector(const ElementType & type, const GhostType & ghost_type, const Vector & vect); inline void free(); }; /// to store data Vector by element type typedef ByElementTypeVector ByElementTypeReal; /// to store data Vector by element type typedef ByElementTypeVector ByElementTypeInt; /// to store data Vector by element type typedef ByElementTypeVector ByElementTypeUInt; /// Map of data of type UInt stored in a mesh typedef std::map *> UIntDataMap; typedef ByElementType ByElementTypeUIntDataMap; /* -------------------------------------------------------------------------- */ /* Mesh */ /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes * Vector, and the connectivity. The connectivity are stored in by element * types. * * To know all the element types present in a mesh you can get the * Mesh::ConnectivityTypeList * * In order to loop on all element you have to loop on all types like this : * @code Mesh::type_iterator it = mesh.firstType(dim, ghost_type); Mesh::type_iterator end = mesh.lastType(dim, ghost_type); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); const Vector & conn = mesh.getConnectivity(*it); for(UInt e = 0; e < nb_element; ++e) { ... } } @endcode */ class Mesh : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const ID & id = "mesh", const MemoryID & memory_id = 0); /// constructor that use an existing nodes coordinates array, by knowing its ID Mesh(UInt spatial_dimension, const ID & nodes_id, const ID & id = "mesh", const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, Vector & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); virtual ~Mesh(); /// @typedef ConnectivityTypeList list of the types present in a Mesh typedef std::set ConnectivityTypeList; // typedef Vector * NormalsMap[_max_element_type]; private: /// initialize the connectivity to NULL and other stuff void init(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /// init a by-element-type real vector with provided ids template void initByElementTypeVector(ByElementTypeVector & v, - UInt nb_component, UInt size, - const bool & flag_nb_node_per_elem_multiply=false, + UInt nb_component, + UInt size, + const bool & flag_nb_node_per_elem_multiply = false, ElementKind element_kind = _ek_regular) const; /// @todo: think about nicer way to do it /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Vector & nodal_values, T * elemental_values, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const; /// extract coordinates of nodes from a reversed element inline void extractNodalCoordinatesFromPBCElement(Real * local_coords, UInt * connectivity, UInt n_nodes); /// convert a element to a linearized element inline UInt elementToLinearized(const Element & elem); /// convert a linearized element to an element inline Element linearizedToElement (UInt linearized_element); /// update the types offsets array for the conversions inline void updateTypesOffsets(const GhostType & ghost_type); /// add a Vector of connectivity for the type . inline void addConnectivityType(const ElementType & type); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const ID &); /// get the spatial dimension of the mesh = number of component of the coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Vector aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Vector &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Vector &); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->getSize(), UInt); /// get the Vector of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Vector &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Vector AKANTU_GET_MACRO(NodesType, *nodes_type, const Vector &); inline Int getNodeType(UInt local_id) const; /// say if a node is a pure ghost node inline bool isPureGhostNode(UInt n) const; /// say if a node is pur local or master node inline bool isLocalOrMasterNode(UInt n) const; inline bool isLocalNode(UInt n) const; inline bool isMasterNode(UInt n) const; inline bool isSlaveNode(UInt n) const; AKANTU_GET_MACRO(XMin, lower_bounds[0], Real); AKANTU_GET_MACRO(YMin, lower_bounds[1], Real); AKANTU_GET_MACRO(ZMin, lower_bounds[2], Real); AKANTU_GET_MACRO(XMax, upper_bounds[0], Real); AKANTU_GET_MACRO(YMax, upper_bounds[1], Real); AKANTU_GET_MACRO(ZMax, upper_bounds[2], Real); inline void getLowerBounds(Real * lower) const; inline void getUpperBounds(Real * upper) const; inline void getLocalLowerBounds(Real * lower) const; inline void getLocalUpperBounds(Real * upper) const; /// get the number of surfaces AKANTU_GET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the connectivity Vector for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); /// @todo take out this set, if mesh can read surface id /// set the number of surfaces AKANTU_SET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; // /// get the number of ghost element of a type in the mesh // inline UInt getNbGhostElement(const ElementType & type) const; /// get the connectivity list either for the elements or the ghost elements inline const ConnectivityTypeList & getConnectivityTypeList(const GhostType & ghost_type = _not_ghost) const; /// compute the barycenter of a given element inline void getBarycenter(UInt element, const ElementType & type, Real * barycenter, GhostType ghost_type = _not_ghost) const; /// get the element connected to a subelement AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementToSubelement, element_to_subelement, Vector); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementToSubelement, element_to_subelement, Vector); /// get the subelement connected to an element AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(SubelementToElement, subelement_to_element, Element); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(SubelementToElement, subelement_to_element, Element); /// get the surface values of facets AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(SurfaceID, surface_id, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(SurfaceID, surface_id, UInt); /// set the int data to the surface id vectors void setSurfaceIDsFromIntData(const std::string & data_name); inline const Vector & getUIntData(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get the kind of the element type static inline ElementKind getKind(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get local connectivity of a facet for a given facet type static inline UInt ** getFacetLocalConnectivity(const ElementType & type); /// get the type of the surface element associated to a given element static inline ElementType getFacetElementType(const ElementType & type); /// get the pointer to the list of elements for a given type inline Vector * getReversedElementsPBCPointer(const ElementType & type); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ typedef ByElementTypeUInt::type_iterator type_iterator; inline type_iterator firstType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.firstType(dim, ghost_type, kind); } inline type_iterator lastType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.lastType(dim, ghost_type, kind); } /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshIOMSH; friend class MeshIOMSHStruct; friend class MeshIODiana; friend class MeshUtils; friend class DistributedSynchronizer; AKANTU_GET_MACRO(NodesPointer, nodes, Vector *); /// get a pointer to the nodes_global_ids Vector and create it if necessary inline Vector * getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Vector and create it if necessary inline Vector * getNodesTypePointer(); /// get a pointer to the connectivity Vector for the given type and create it if necessary inline Vector * getConnectivityPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); // inline Vector * getNormalsPointer(ElementType type) const; /// get a pointer to the surface_id Vector for the given type and create it if necessary inline Vector * getSurfaceIDPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get the UIntDataMap for a given ElementType inline UIntDataMap & getUIntDataMap(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the IntDataMap pointer (modifyable) for a given ElementType inline Vector * getUIntDataPointer(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type = _not_ghost); /// get a pointer to the element_to_subelement Vector for the given type and create it if necessary inline Vector > * getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Vector for the given type and create it if necessary inline Vector * getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// id of the mesh ID id; /// array of the nodes coordinates Vector * nodes; /// global node ids Vector * nodes_global_ids; /// node type, -3 pure ghost, -2 master for the node, -1 normal node, i in /// [0-N] slave node and master is proc i Vector * nodes_type; /// global number of nodes; UInt nb_global_nodes; /// boolean to know if the nodes have to be deleted with the mesh or not bool created_nodes; /// all class of elements present in this mesh (for heterogenous meshes) ByElementTypeUInt connectivities; /// map to normals for all class of elements present in this mesh ByElementTypeReal normals; /// list of all existing types in the mesh ConnectivityTypeList type_set; /// the spatial dimension of this mesh UInt spatial_dimension; /// types offsets Vector types_offsets; /// list of all existing types in the mesh ConnectivityTypeList ghost_type_set; /// ghost types offsets Vector ghost_types_offsets; /// number of surfaces present in this mesh UInt nb_surfaces; /// surface id of the surface elements in this mesh ByElementTypeUInt surface_id; /// min of coordinates Real lower_bounds[3]; /// max of coordinates Real upper_bounds[3]; /// size covered by the mesh on each direction Real size[3]; /// local min of coordinates Real local_lower_bounds[3]; /// local max of coordinates Real local_upper_bounds[3]; /// List of elements connected to subelements ByElementTypeVector > element_to_subelement; /// List of subelements connected to elements ByElementTypeVector subelement_to_element; // /// list of elements that are reversed due to pbc // ByElementTypeUInt reversed_elements_pbc; // /// direction in which pbc are to be applied // UInt pbc_directions[3]; /// list of the vectors corresponding to tags in the mesh ByElementTypeUIntDataMap uint_data; }; /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Element & _this) { _this.printself(stream); return stream; } #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "mesh_inline_impl.cc" #endif #include "by_element_type_tmpl.hh" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MESH_HH__ */ diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index 60392a4bd..e391d9ed4 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,780 +1,795 @@ /** * @file material.cc * @author Nicolas Richart * @date Tue Jul 27 11:43:41 2010 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(model.getMemoryID()), id(id), name(""), model(&model), stress("stress", id), strain("strain", id), element_filter("element_filter", id), // potential_energy_vector(false), potential_energy("potential_energy", id), is_non_local(false), - interpolationInvCoordinates("interpolationInvCoordinates", id), + interpolation_inverse_coordinates("interpolation inverse coordinates", id), + interpolation_points_matrices("interpolation points matrices", id), is_init(false) { AKANTU_DEBUG_IN(); rho = 0; AKANTU_DEBUG_ASSERT(this->model,"model has wrong type: cannot proceed"); spatial_dimension = this->model->getSpatialDimension(); /// allocate strain stress for local elements initInternalVector(strain, spatial_dimension * spatial_dimension); initInternalVector(stress, spatial_dimension * spatial_dimension); /// for each connectivity types allocate the element filer array of the material initInternalVector(element_filter, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool Material::setParam(const std::string & key, const std::string & value, __attribute__ ((unused)) const ID & id) { std::stringstream sstr(value); if(key == "name") name = std::string(value); else if(key == "rho") { sstr >> rho; } else return false; return true; } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); resizeInternalVector(stress); resizeInternalVector(strain); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const { AKANTU_DEBUG_IN(); model->getFEM().getMesh().initByElementTypeVector(vect, nb_component, spatial_dimension, false, element_kind); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::resizeInternalVector(ByElementTypeVector & by_el_type_vect, ElementKind element_kind) const { AKANTU_DEBUG_IN(); FEM * fem = & model->getFEM(); if (element_kind == _ek_cohesive) fem = & model->getFEM("CohesiveFEM"); const Mesh & mesh = fem->getMesh(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, element_kind); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, element_kind); for(; it != end; ++it) { const Vector & elem_filter = element_filter(*it, gt); UInt nb_element = elem_filter.getSize(); UInt nb_quadrature_points = fem->getNbQuadraturePoints(*it, gt); UInt new_size = nb_element * nb_quadrature_points; Vector & vect = by_el_type_vect(*it, gt); UInt size = vect.getSize(); UInt nb_component = vect.getNbComponent(); vect.resize(new_size); memset(vect.storage() + size * nb_component, 0, (new_size - size) * nb_component * sizeof(T)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::updateResidual(Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); computeAllStresses(displacement, ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Vector & residual = const_cast &>(model->getResidual()); Mesh & mesh = model->getFEM().getMesh(); 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) { const Vector & shapes_derivatives = model->getFEM().getShapesDerivatives(*it, ghost_type); Vector & elem_filter = element_filter(*it, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, ghost_type); UInt nb_element = elem_filter.getSize(); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Vector * sigma_dphi_dx = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Real * shapesd = shapes_derivatives.storage(); Real * shapesd_val; UInt * elem_filter_val = elem_filter.storage(); Vector * shapesd_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "filtered shapesd"); Real * shapesd_filtered_val = shapesd_filtered->values; for (UInt el = 0; el < nb_element; ++el) { shapesd_val = shapesd + elem_filter_val[el] * size_of_shapes_derivatives * nb_quadrature_points; memcpy(shapesd_filtered_val, shapesd_val, size_of_shapes_derivatives * nb_quadrature_points * sizeof(Real)); shapesd_filtered_val += size_of_shapes_derivatives * nb_quadrature_points; } Vector & stress_vect = stress(*it, ghost_type); // Vector::iterator sigma = stress_vect.begin(spatial_dimension, spatial_dimension); // Vector::iterator sigma_end = stress_vect.end(spatial_dimension, spatial_dimension); // Vector::iterator nabla_B = shapesd_filtered->begin(nb_nodes_per_element, spatial_dimension); // Vector::iterator sigma_dphi_dx_it = sigma_dphi_dx->begin(nb_nodes_per_element, spatial_dimension); // for (; sigma != sigma_end; ++sigma, ++nabla_B, ++sigma_dphi_dx_it) { // sigma_dphi_dx_it->mul(*nabla_B, *sigma); // } Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, *shapesd_filtered, stress_vect, *sigma_dphi_dx); delete shapesd_filtered; /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ Vector * int_sigma_dphi_dx = new Vector(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); model->getFEM().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, *it, ghost_type, &elem_filter); delete sigma_dphi_dx; /// assemble model->getFEM().assembleVector(*int_sigma_dphi_dx, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, &elem_filter, -1); delete int_sigma_dphi_dx; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the strain * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); resizeInternalVector(stress); resizeInternalVector(strain); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); Vector & strain_vect = strain(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(displacement, strain_vect, spatial_dimension, *it, ghost_type, &elem_filter); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Vector & displacement = model->getDisplacement(); resizeInternalVector(strain); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); Vector & strain_vect = strain(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(displacement, strain_vect, spatial_dimension, *it, ghost_type, &elem_filter); setToSteadyState(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(Vector & current_position, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh & mesh = model->getFEM().getMesh(); 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) { switch(spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(current_position, *it, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(current_position, *it, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(current_position, *it, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(Vector & current_position, const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); const Vector & shapes_derivatives = model->getFEM().getShapesDerivatives(type,ghost_type); Vector & elem_filter = element_filter(type, ghost_type); Vector & strain_vect = strain(type, ghost_type); UInt * elem_filter_val = elem_filter.storage(); UInt nb_element = elem_filter.getSize(); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type); strain_vect.resize(nb_quadrature_points * nb_element); model->getFEM().gradientOnQuadraturePoints(current_position, strain_vect, dim, type, ghost_type, &elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Vector * tangent_stiffness_matrix = new Vector(nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); computeTangentStiffness(type, *tangent_stiffness_matrix, ghost_type); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Vector * bt_d_b = new Vector(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); UInt size_of_b = tangent_size * bt_d_b_size; Real * B = new Real[size_of_b]; Real * Bt_D = new Real[size_of_b]; Real * Bt_D_B = bt_d_b->storage(); Real * D = tangent_stiffness_matrix->storage(); UInt offset_bt_d_b = bt_d_b_size * bt_d_b_size; UInt offset_d = tangent_size * tangent_size; for (UInt e = 0; e < nb_element; ++e) { Real * shapes_derivatives_val = shapes_derivatives.values + elem_filter_val[e]*size_of_shapes_derivatives*nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q) { transferBMatrixToSymVoigtBMatrix(shapes_derivatives_val, B, nb_nodes_per_element); Math::matrixt_matrix(bt_d_b_size, tangent_size, tangent_size, B, D, Bt_D); Math::matrix_matrix(bt_d_b_size, bt_d_b_size, tangent_size, Bt_D, B, Bt_D_B); shapes_derivatives_val += size_of_shapes_derivatives; D += offset_d; Bt_D_B += offset_bt_d_b; } } delete [] B; delete [] Bt_D; delete tangent_stiffness_matrix; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Vector * K_e = new Vector(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); model->getFEM().integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, &elem_filter); delete bt_d_b; model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, &elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(ghost_type != _not_ghost) return; Real * epot = potential_energy(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; computePotentialEnergyOnQuad(grad_u, sigma, *epot); epot++; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElement() { AKANTU_DEBUG_IN(); Mesh & mesh = model->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(!potential_energy.exists(*it, _not_ghost)) { UInt nb_element = element_filter(*it, _not_ghost).getSize(); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, _not_ghost); potential_energy.alloc(nb_element * nb_quadrature_points, 1, *it, _not_ghost); } computePotentialEnergy(*it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElement(); /// integrate the potential energy for each type of elements Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { epot += model->getFEM().integrate(potential_energy(*it, _not_ghost), *it, _not_ghost, &element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if(type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates) const { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEM().getMesh(); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); Vector nodes_coordinates(mesh.getNodes(), true); nodes_coordinates += model->getDisplacement(); for(UInt gt = (UInt) _not_ghost; gt < (UInt) _casper; ++gt) { GhostType ghost_type = (GhostType) gt; 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) { const Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_tot_quad = model->getFEM().getNbQuadraturePoints(*it, ghost_type) * nb_element; Vector & quads = quadrature_points_coordinates(*it, ghost_type); quads.resize(nb_tot_quad); model->getFEM().interpolateOnQuadraturePoints(nodes_coordinates, quads, spatial_dimension, *it, ghost_type, &elem_filter); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void Material::initInterpolateElementalField() { +void Material::initElementalFieldInterpolation(ByElementTypeReal & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); - const Mesh & mesh = model->getFEM().getMesh(); - const Vector & position = mesh.getNodes(); + ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", id); + computeQuadraturePointsCoordinates(quadrature_points_coordinates); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); - for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element == 0) continue; - - /// compute quadrature point position - UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(*it); - UInt nb_quad = nb_quad_per_element * nb_element; - Vector quad_coordinates(nb_quad, spatial_dimension); - - model->getFEM().interpolateOnQuadraturePoints(position, - quad_coordinates, - spatial_dimension, - *it); - - interpolationInvCoordinates.alloc(nb_quad, nb_quad_per_element, *it); - - Vector & interp_inv_coord = interpolationInvCoordinates(*it); - ElementType type = *it; +#define AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD(type) \ + initElementalFieldInterpolation(quadrature_points_coordinates(type), \ + interpolation_points_coordinates(type)) -#define INIT_INTERPOLATE_ELEMENTAL_FIELD(type) \ - initInterpolation(quad_coordinates, interp_inv_coord) \ - - AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_INTERPOLATE_ELEMENTAL_FIELD); -#undef INIT_INTERPOLATE_ELEMENTAL_FIELD + AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD); +#undef AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template -void Material::initInterpolation(const Vector & quad_coordinates, - Vector & interp_inv_coord) { +void Material::initElementalFieldInterpolation(const Vector & quad_coordinates, + const Vector & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); + UInt size_inverse_coords = getSizeElementalFieldInterpolationCoodinates(); Vector & elem_fil = element_filter(type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type); - types::Matrix quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); + UInt nb_interpolation_points = interpolation_points_coordinates.getSize(); + AKANTU_DEBUG_ASSERT(nb_interpolation_points % nb_element == 0, + "Can't interpolate elemental field on elements, the coordinates vector has a wrong size"); + UInt nb_interpolation_points_per_elem = nb_interpolation_points / nb_element; + + if(!interpolation_inverse_coordinates.exists(type)) + interpolation_inverse_coordinates.alloc(nb_element, + size_inverse_coords*size_inverse_coords, + type); + + if(!interpolation_points_matrices.exists(type)) + interpolation_points_matrices.alloc(nb_element, + nb_interpolation_points_per_elem * size_inverse_coords, + type); + + Vector & interp_inv_coord = interpolation_inverse_coordinates(type); + Vector & interp_points_mat = interpolation_points_matrices(type); + + types::Matrix quad_coord_matrix(size_inverse_coords, size_inverse_coords); - UInt quad_block = nb_quad_per_element * spatial_dimension; - UInt nb_quad2 = nb_quad_per_element * nb_quad_per_element; + Vector::const_iterator quad_coords_it = + quad_coordinates.begin_reinterpret(nb_quad_per_element, + spatial_dimension, + nb_element, + nb_quad_per_element * spatial_dimension); + Vector::const_iterator points_coords_it = + interpolation_points_coordinates.begin_reinterpret(nb_interpolation_points_per_elem, + spatial_dimension, + nb_element, + nb_interpolation_points_per_elem * spatial_dimension); + + Vector::iterator inv_quad_coord_it = + interp_inv_coord.begin(size_inverse_coords, size_inverse_coords); + + Vector::iterator inv_points_mat_it = + interp_points_mat.begin(nb_interpolation_points_per_elem, size_inverse_coords); + /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el) { + /// matrix containing the quadrature points coordinates + const types::Matrix & quad_coords = *quad_coords_it; + /// matrix to store the matrix inversion result + types::Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; + + /// insert the quad coordinates in a matrix compatible with the interpolation + buildElementalFieldInterpolationCoodinates(quad_coords, + quad_coord_matrix); + + /// invert the interpolation matrix + inv_quad_coord_matrix.inverse(quad_coord_matrix); + - /// create a matrix containing the quadrature points coordinates - types::Matrix quad_coords(quad_coordinates.storage() - + quad_block * elem_fil(el), - nb_quad_per_element, - spatial_dimension); + /// matrix containing the interpolation points coordinates + const types::Matrix & points_coords = *points_coords_it; + /// matrix to store the interpolation points coordinates + /// compatible with these functions + types::Matrix & inv_points_coord_matrix = *inv_points_mat_it; /// insert the quad coordinates in a matrix compatible with the interpolation - buildInterpolationCoodinates(quad_coords, quad_coord_matrix); + buildElementalFieldInterpolationCoodinates(points_coords, + inv_points_coord_matrix); - /// create a matrix to store the matrix inversion result - types::Matrix inv_quad_coord_matrix(interp_inv_coord.storage() - + nb_quad2 * el, - nb_quad_per_element, - nb_quad_per_element); - /// invert the interpolation matrix - Math::inv(nb_quad_per_element, - quad_coord_matrix.storage(), - inv_quad_coord_matrix.storage()); + ++inv_quad_coord_it; + ++inv_points_mat_it; + ++quad_coords_it; + ++points_coords_it; } - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(const ElementType type, - const Vector & coordinates, Vector & result) { AKANTU_DEBUG_IN(); #define INTERPOLATE_ELEMENTAL_FIELD(type) \ interpolateElementalField(stress(type), \ - coordinates, \ result) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTERPOLATE_ELEMENTAL_FIELD); #undef INTERPOLATE_ELEMENTAL_FIELD AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::interpolateElementalField(const Vector & field, - const Vector & coordinates, Vector & result) { AKANTU_DEBUG_IN(); Vector & elem_fil = element_filter(type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type); + UInt size_inverse_coords = getSizeElementalFieldInterpolationCoodinates(); types::Matrix coefficients(nb_quad_per_element, field.getNbComponent()); - const Vector & interp_inv_coord = interpolationInvCoordinates(type); + const Vector & interp_inv_coord = interpolation_inverse_coordinates(type); + const Vector & interp_points_coord = interpolation_points_matrices(type); - Vector::const_iterator field_it = - field.begin_reinterpret(nb_quad_per_element, - field.getNbComponent(), - nb_element, - nb_quad_per_element * field.getNbComponent()); + UInt nb_interpolation_points_per_elem = interp_points_coord.getNbComponent() / size_inverse_coords; - AKANTU_DEBUG_ASSERT(coordinates.getSize() % nb_element == 0, - "Can't interpolate elemental field on elements, the coordinates vector has a wrong size"); + Vector filtered_field(nb_element, nb_quad_per_element); + extractElementalFieldForInterplation(field, filtered_field); - UInt nb_points_per_elem = coordinates.getSize() / nb_element; + Vector::const_iterator field_it + = field.begin_reinterpret(nb_quad_per_element, + field.getNbComponent(), + nb_element, + nb_quad_per_element * field.getNbComponent()); - types::Matrix coord_matrix(nb_points_per_elem, nb_quad_per_element); + Vector::const_iterator interpolation_points_coordinates_it = + interp_points_coord.begin(nb_interpolation_points_per_elem, size_inverse_coords); Vector::iterator result_it - = result.begin_reinterpret(nb_points_per_elem, + = result.begin_reinterpret(nb_interpolation_points_per_elem, field.getNbComponent(), nb_element, - nb_points_per_elem * field.getNbComponent()); + nb_interpolation_points_per_elem * field.getNbComponent()); - UInt coord_block = nb_points_per_elem * spatial_dimension; - UInt nb_quad2 = nb_quad_per_element * nb_quad_per_element; + Vector::const_iterator inv_quad_coord_it = + interp_inv_coord.begin(size_inverse_coords, size_inverse_coords); /// loop over the elements of the current material and element type - for (UInt el = 0; el < nb_element; ++el, ++field_it, ++result_it) { - + for (UInt el = 0; el < nb_element; + ++el, ++field_it, ++result_it, + ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { /** - * create a matrix containing the inversion of the quadrature - * points' coordinates - * + * matrix containing the inversion of the quadrature points' + * coordinates */ - types::Matrix inv_quad_coord_matrix(interp_inv_coord.storage() - + nb_quad2 * el, - nb_quad_per_element, - nb_quad_per_element); - + const types::Matrix & inv_quad_coord_matrix = *inv_quad_coord_it; /** * multiply it by the field values over quadrature points to get * the interpolation coefficients - * */ coefficients.mul(inv_quad_coord_matrix, *field_it); - /// create a matrix containing the points' coordinates - types::Matrix coord(coordinates.storage() - + coord_block * elem_fil(el), - nb_points_per_elem, - spatial_dimension); - - /// insert the points coordinates in a matrix compatible with the interpolation - buildInterpolationCoodinates(coord, coord_matrix); + /// matrix containing the points' coordinates + const types::Matrix & coord = *interpolation_points_coordinates_it; /// multiply the coordinates matrix by the coefficients matrix and store the result - (*result_it).mul(coord_matrix, coefficients); + (*result_it).mul(coord, coefficients); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const Vector & Material::getVector(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << id << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getVector(fvect_id); } catch(debug::Exception & e) { AKANTU_EXCEPTION("The material " << name << "(" <(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const; template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const; template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; __END_AKANTU__ diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 38e1256ff..ea49175ff 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,495 +1,503 @@ /** * @file material.hh * @author Nicolas Richart * @date Fri Jul 23 09:06:29 2010 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" //#include "fem.hh" //#include "mesh.hh" #include "data_accessor.hh" //#include "static_communicator.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; class CommunicationBuffer; } __BEGIN_AKANTU__ /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(const ElementType & el_type, * Vector & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : protected Memory, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const ID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(Vector & displacement, GhostType ghost_type = _not_ghost); void assembleResidual(GhostType ghost_type); /// compute the residual for this material virtual void computeAllStresses(Vector & current_position, GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(Vector & current_position, GhostType ghost_type); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(Real h, const Element & element = ElementNull) = 0; /// compute the p-wave speed in the material virtual Real getPushWaveSpeed() { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// compute the s-wave speed in the material virtual Real getShearWaveSpeed() { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; - /// interpolate stress on given positions for each element by means - /// of a geometrical interpolation on quadrature points + /** + * interpolate stress on given positions for each element by means + * of a geometrical interpolation on quadrature points + */ virtual void interpolateStress(const ElementType type, - const Vector & coordinates, Vector & result); - /// function to initialize the elemental field interpolation - /// function by inverting the quadrature points' coordinates - virtual void initInterpolateElementalField(); + /** + * function to initialize the elemental field interpolation + * function by inverting the quadrature points' coordinates + */ + virtual void initElementalFieldInterpolation(ByElementTypeReal & interpolation_points_coordinates); protected: /// constitutive law virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// set the material to steady state (to be implemented for materials that need it) virtual void setToSteadyState(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {}; - // /// constitutive law - // virtual void computeNonLocalStress(ElementType el_type, - // GhostType ghost_type = _not_ghost) { - // AKANTU_DEBUG_TO_IMPLEMENT(); - // }; - /// compute the tangent stiffness matrix virtual void computeTangentStiffness(const ElementType & el_type, Vector & tangent_matrix, GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); template void assembleStiffnessMatrix(Vector & current_position, const ElementType & type, GhostType ghost_type); /// transfer the B matrix to a Voigt notation B matrix template inline void transferBMatrixToSymVoigtBMatrix(Real * B, Real * Bvoigt, UInt nb_nodes_per_element) const; inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// compute the potential energy by element void computePotentialEnergyByElement(); - + /// compute the coordinates of the quadrature points void computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates) const; /// interpolate an elemental field on given points for each element template void interpolateElementalField(const Vector & field, - const Vector & coordinates, Vector & result); /// template function to initialize the elemental field interpolation template - void initInterpolation(const Vector & quad_coordinates, - Vector & interp_inv_coord); + void initElementalFieldInterpolation(const Vector & quad_coordinates, + const Vector & interpolation_points_coordinates); /// build the coordinate matrix for the interpolation on elemental field template - inline void buildInterpolationCoodinates(const types::Matrix & coordinates, - types::Matrix & coordMatrix); + inline void buildElementalFieldInterpolationCoodinates(const types::Matrix & coordinates, + types::Matrix & coordMatrix); + + /// get the size of the coordiante matrix used in the interpolation + template + inline UInt getSizeElementalFieldInterpolationCoodinates(); + + /// extract the field values corresponding to the quadrature points used for the interpolation + template + inline void extractElementalFieldForInterplation(const Vector & field, + Vector & filtered_field); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// compute the potential energy for on element inline void computePotentialEnergyOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & epot); - + public: /// allocate an internal vector template void initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind = _ek_regular) const; /// resize an internal vector template void resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind = _ek_regular) const; /* ------------------------------------------------------------------------ */ template inline void gradUToF(const types::Matrix & grad_u, types::Matrix & F); inline void rightCauchy(const types::Matrix & F, types::Matrix & C); inline void leftCauchy (const types::Matrix & F, types::Matrix & B); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataToPack(__attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline UInt getNbDataToUnpack(__attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual UInt getNbDataToPack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual UInt getNbDataToUnpack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline void packData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { } virtual void packData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const UInt index, __attribute__((unused)) SynchronizationTag tag) const { } virtual inline void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) { } virtual void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const UInt index, __attribute__((unused)) SynchronizationTag tag) { } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); Real getPotentialEnergy(); virtual Real getEnergy(std::string type); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Strain, strain, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); const Vector & getVector(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; virtual Real getProperty(const ID & param) const; virtual void setProperty(const ID & param, Real value); protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material ID id; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// stresses arrays ordered by element types ByElementTypeReal stress; /// strains arrays ordered by element types ByElementTypeReal strain; /// list of element handled by the material ByElementTypeUInt element_filter; /// is the vector for potential energy initialized // bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// tell if using in non local mode or not bool is_non_local; /// spatial dimension UInt spatial_dimension; - /// elemental field interpolation coefficients - ByElementTypeReal interpolationInvCoordinates; + /// elemental field interpolation coordinates + ByElementTypeReal interpolation_inverse_coordinates; + + /// elemental field interpolation points + ByElementTypeReal interpolation_points_matrices; private: /// boolean to know if the material has been initialized bool is_init; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "material_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ Vector::iterator stress_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ for(;strain_it != strain_end; ++strain_it, ++stress_it) { \ types::Matrix & __attribute__((unused)) grad_u = *strain_it; \ types::Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Vector::iterator tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;strain_it != strain_end; ++strain_it, ++tangent_it) { \ types::Matrix & __attribute__((unused)) grad_u = *strain_it; \ types::Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } /* -------------------------------------------------------------------------- */ /* Material list */ /* -------------------------------------------------------------------------- */ #define AKANTU_CORE_MATERIAL_LIST \ ((2, (elastic , MaterialElastic ))) \ ((2, (viscoelastic , MaterialViscoElastic ))) \ ((2, (elastic_orthotropic , MaterialElasticOrthotropic ))) \ ((2, (elastic_caughey , MaterialElasticCaughey ))) \ ((2, (neohookean , MaterialNeohookean ))) #define AKANTU_COHESIVE_MATERIAL_LIST \ ((2, (cohesive_bilinear , MaterialCohesiveBilinear ))) \ ((2, (cohesive_linear , MaterialCohesiveLinear ))) \ ((2, (cohesive_linear_extrinsic, MaterialCohesiveLinearExtrinsic ))) \ ((2, (cohesive_linear_exponential_extrinsic, MaterialCohesiveLinearExponentialExtrinsic ))) \ ((2, (cohesive_exponential , MaterialCohesiveExponential ))) #define AKANTU_DAMAGE_MATERIAL_LIST \ ((2, (damage_linear , MaterialDamageLinear ))) \ ((2, (marigo , MaterialMarigo ))) \ ((2, (mazars , MaterialMazars ))) \ ((2, (vreepeerlings , MaterialVreePeerlings ))) #ifdef AKANTU_DAMAGE_NON_LOCAL #define AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST \ ((stress_wf, StressBasedWeightFunction )) \ ((damage_wf, DamagedWeightFunction )) \ ((remove_wf, RemoveDamagedWeightFunction)) \ ((base_wf, BaseWeightFunction )) # define AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST \ ((3, (marigo_non_local , MaterialMarigoNonLocal, AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST))) \ ((2, (mazars_non_local , MaterialMazarsNonLocal ))) \ ((3, (vreepeerlings_non_local, MaterialVreePeerlingsNonLocal, AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST))) #else # define AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST #endif #define AKANTU_MATERIAL_LIST \ AKANTU_CORE_MATERIAL_LIST \ AKANTU_COHESIVE_MATERIAL_LIST \ AKANTU_DAMAGE_MATERIAL_LIST \ AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST #define INSTANSIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #if defined(__INTEL_COMPILER) #pragma warning ( push ) /* warning #654: overloaded virtual function "akantu::Material::computeStress" is only partially overridden in class "akantu::Material*" */ #pragma warning ( disable : 654 ) #endif //defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ // elastic materials #include "material_elastic.hh" #include "material_elastic_caughey.hh" #include "material_viscoelastic.hh" #include "material_neohookean.hh" #include "material_elastic_orthotropic.hh" // damage materials #include "material_damage.hh" #include "material_marigo.hh" #include "material_mazars.hh" #include "material_damage_linear.hh" #include "material_vreepeerlings.hh" #if defined(AKANTU_DAMAGE_NON_LOCAL) # include "material_marigo_non_local.hh" # include "material_mazars_non_local.hh" # include "material_vreepeerlings_non_local.hh" #endif // cohesive materials #include "material_cohesive.hh" #include "material_cohesive_linear.hh" #include "material_cohesive_bilinear.hh" #include "material_cohesive_linear_extrinsic.hh" #include "material_cohesive_exponential.hh" #include "material_cohesive_linear_exponential_extrinsic.hh" #if defined(__INTEL_COMPILER) #pragma warning ( pop ) #endif //defined(__INTEL_COMPILER) #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/material_inline_impl.cc b/src/model/solid_mechanics/material_inline_impl.cc index db7c02212..20e3e5aaf 100644 --- a/src/model/solid_mechanics/material_inline_impl.cc +++ b/src/model/solid_mechanics/material_inline_impl.cc @@ -1,160 +1,180 @@ /** * @file material_inline_impl.cc * @author Nicolas Richart * @date Tue Jul 27 11:57:43 2010 * * @brief Implementation of the inline functions of the class material * * @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 . * */ /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const ElementType & type, UInt element, const GhostType & ghost_type) { element_filter(type, ghost_type).push_back(element); return element_filter(type, ghost_type).getSize()-1; } /* -------------------------------------------------------------------------- */ inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const { return (dim * (dim - 1) / 2 + dim); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToF(const types::Matrix & grad_u, types::Matrix & F) { UInt size_F = F.size(); AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim, "The dimension of the tensor F should be greater or equal to the dimension of the tensor grad_u."); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) F(i, j) = grad_u(i, j); for (UInt i = 0; i < size_F; ++i) F(i, i) += 1; } /* -------------------------------------------------------------------------- */ inline void Material::rightCauchy(const types::Matrix & F, types::Matrix & C) { C.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::leftCauchy(const types::Matrix & F, types::Matrix & B) { B.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::computePotentialEnergyOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & epot) { epot = 0.; for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) epot += sigma(i, j) * grad_u(i, j); epot *= .5; } /* -------------------------------------------------------------------------- */ template inline void Material::transferBMatrixToSymVoigtBMatrix(Real * B, Real * Bvoigt, UInt nb_nodes_per_element) const { UInt size = getTangentStiffnessVoigtSize(dim) * nb_nodes_per_element * dim; memset(Bvoigt, 0, size * sizeof(Real)); for (UInt i = 0; i < dim; ++i) { Real * Bvoigt_tmp = Bvoigt + i * (dim * nb_nodes_per_element + 1); Real * Btmp = B + i; for (UInt n = 0; n < nb_nodes_per_element; ++n) { *Bvoigt_tmp = *Btmp; Btmp += dim; Bvoigt_tmp += dim; } } if(dim == 2) { ///in 2D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{\partial N_i}{\partial y}]@f$ row Real * Bvoigt_tmp = Bvoigt + dim * nb_nodes_per_element * 2; for (UInt n = 0; n < nb_nodes_per_element; ++n) { Bvoigt_tmp[1] = B[n * dim + 0]; Bvoigt_tmp[0] = B[n * dim + 1]; Bvoigt_tmp += dim; } } if(dim == 3) { UInt Bvoigt_wcol = dim * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { Real dndx = B[n * dim + 0]; Real dndy = B[n * dim + 1]; Real dndz = B[n * dim + 2]; UInt Bni_off = n * dim; ///in 3D, fill the @f$ [0, \frac{\partial N_i}{\partial y}, \frac{N_i}{\partial z}]@f$ row Bvoigt[3 * Bvoigt_wcol + Bni_off + 1] = dndz; Bvoigt[3 * Bvoigt_wcol + Bni_off + 2] = dndy; ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, 0, \frac{N_i}{\partial z}]@f$ row Bvoigt[4 * Bvoigt_wcol + Bni_off + 0] = dndz; Bvoigt[4 * Bvoigt_wcol + Bni_off + 2] = dndx; ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{N_i}{\partial y}, 0]@f$ row Bvoigt[5 * Bvoigt_wcol + Bni_off + 0] = dndy; Bvoigt[5 * Bvoigt_wcol + Bni_off + 1] = dndx; } } } /* -------------------------------------------------------------------------- */ template -inline void Material::buildInterpolationCoodinates(__attribute__((unused)) const types::Matrix & coordinates, - __attribute__((unused)) types::Matrix & coordMatrix) { +inline void Material::buildElementalFieldInterpolationCoodinates(__attribute__((unused)) const types::Matrix & coordinates, + __attribute__((unused)) types::Matrix & coordMatrix) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template<> -inline void Material::buildInterpolationCoodinates<_triangle_3>(const types::Matrix & coordinates, - types::Matrix & coordMatrix) { +inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_3>(const types::Matrix & coordinates, + types::Matrix & coordMatrix) { for (UInt i = 0; i < coordinates.rows(); ++i) coordMatrix(i, 0) = 1; } /* -------------------------------------------------------------------------- */ template<> -inline void Material::buildInterpolationCoodinates<_triangle_6>(const types::Matrix & coordinates, - types::Matrix & coordMatrix) { +inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_6>(const types::Matrix & coordinates, + types::Matrix & coordMatrix) { UInt nb_quadrature_points = 3; for (UInt i = 0; i < coordinates.rows(); ++i) { coordMatrix(i, 0) = 1; for (UInt j = 1; j < nb_quadrature_points; ++j) coordMatrix(i, j) = coordinates(i, j-1); } } + +__END_AKANTU__ + +#include "solid_mechanics_model.hh" + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +template +inline UInt Material::getSizeElementalFieldInterpolationCoodinates() { + return model->getFEM().getNbQuadraturePoints(type); +} + +/* -------------------------------------------------------------------------- */ +template +inline void Material::extractElementalFieldForInterplation(const Vector & field, + Vector & filtered_field) { + filtered_field.copy(field); +} + diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc index 2f8b8d450..b3c41b08f 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc @@ -1,344 +1,342 @@ /** * @file material_cohesive_linear.cc * @author Marco Vocialta * @date Mon Feb 20 12:14:16 2012 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear_extrinsic.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearExtrinsic::MaterialCohesiveLinearExtrinsic(SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model,id), sigma_c_eff("sigma_c_eff",id), delta_c("delta_c",id), normal_stress(spatial_dimension), tangential_stress(spatial_dimension) { AKANTU_DEBUG_IN(); sigma_c = 0; beta = 0; G_cI = 0; G_cII = 0; rand = 0; penalty = 0; initInternalVector(sigma_c_eff, 1, _ek_cohesive); initInternalVector(delta_c, 1, _ek_cohesive); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearExtrinsic::~MaterialCohesiveLinearExtrinsic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); kappa = G_cII / G_cI; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::resizeCohesiveVectors() { MaterialCohesive::resizeCohesiveVectors(); resizeInternalVector(sigma_c_eff, _ek_cohesive); resizeInternalVector(delta_c, _ek_cohesive); FEM & fem_cohesive = model->getFEM("CohesiveFEM"); const Mesh & mesh = fem_cohesive.getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { const Vector & elem_filter = element_filter(*it, _not_ghost); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; UInt nb_quadrature_points = fem_cohesive.getNbQuadraturePoints(*it, _not_ghost); UInt nb_element_old = nb_element - sigma_insertion.getSize() / nb_quadrature_points; Vector & sigma_c_eff_vec = sigma_c_eff(*it, _not_ghost); Vector & delta_c_vec = delta_c(*it, _not_ghost); for (UInt el = nb_element_old; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Real new_sigma = sigma_insertion((el - nb_element_old) * nb_quadrature_points+q); Real new_delta = 2 * G_cI / new_sigma; sigma_c_eff_vec(el * nb_quadrature_points + q) = new_sigma; delta_c_vec(el * nb_quadrature_points + q) = new_delta; } } } } /* -------------------------------------------------------------------------- */ template bool MaterialCohesiveLinearExtrinsic::setParam(const std::string & key, const std::string & value, const ID & id) { std::stringstream sstr(value); if(key == "sigma_c") { sstr >> sigma_c; } else if(key == "beta") { sstr >> beta; } else if(key == "G_cI") { sstr >> G_cI; } else if(key == "G_cII") { sstr >> G_cII; } else if(key == "rand") { sstr >> rand; } else if(key == "penalty") { sstr >> penalty; } else { return Material::setParam(key, value, id); } return true; } +/* -------------------------------------------------------------------------- */ +template +inline Real MaterialCohesiveLinearExtrinsic::computeEffectiveNorm(const types::Matrix & stress, + const types::RVector & normal, + const types::RVector & tangent) { + AKANTU_DEBUG_IN(); + + normal_stress.mul(stress, normal); + tangential_stress.mul(stress, tangent); + + Real normal_contrib = normal_stress.dot(normal); + Real tangent_contrib = tangential_stress.dot(tangent); + + normal_contrib = std::max(0., normal_contrib); + + AKANTU_DEBUG_OUT(); + + return std::sqrt(normal_contrib*normal_contrib + + tangent_contrib*tangent_contrib/beta/beta); +} + /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::computeStressNorms(const Vector & facet_stress, - types::RVector & stress_check) { + Vector & stress_check) { AKANTU_DEBUG_IN(); sigma_insertion.resize(0); Vector & facets_check = model->getFacetsCheck(); ElementType type_facet = model->getFacetType(); UInt nb_quad_facet = model->getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_facet = facets_check.getSize(); const Vector & tangents = model->getTangents(); const Vector & normals = model->getFEM("FacetsFEM").getNormalsOnQuadPoints(type_facet); - UInt sp2 = spatial_dimension * spatial_dimension; + Vector::iterator stress_check_it = + stress_check.begin(nb_quad_facet); - for (UInt f = 0; f < nb_facet; ++f) { - if (facets_check(f) == true) { + Vector::const_iterator normal_it = + normals.begin(spatial_dimension); - stress_check.clear(); + Vector::const_iterator tangent_it = + tangents.begin(spatial_dimension); - for (UInt e = 0; e < 2; ++e) { - for (UInt q = 0; q < nb_quad_facet; ++q) { - types::Matrix stress_edge(facet_stress.storage() - + f * 2 * nb_quad_facet * sp2 - + e * nb_quad_facet * sp2 - + q * sp2, - spatial_dimension, spatial_dimension); + Vector::const_iterator facet_stress_it = + facet_stress.begin(spatial_dimension, spatial_dimension); - types::RVector normal(normals.storage() + f*nb_quad_facet*spatial_dimension, - spatial_dimension); + for (UInt f = 0; f < nb_facet; ++f, ++stress_check_it) { + for (UInt q = 0; q < nb_quad_facet; ++q, ++normal_it, ++tangent_it) { + for (UInt e = 0; e < 2; ++e, ++facet_stress_it) { - types::RVector tangent(tangents.storage() + f*nb_quad_facet*spatial_dimension, - spatial_dimension); + if (facets_check(f) == true) { + Real effective_norm = + computeEffectiveNorm(*facet_stress_it, *normal_it, *tangent_it); - stress_check(q) = std::max(stress_check(q), - computeEffectiveNorm(stress_edge, normal, tangent)); + (*stress_check_it)(q) = + std::max((*stress_check_it)(q), effective_norm); } + } } } AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -template -inline Real MaterialCohesiveLinearExtrinsic::computeEffectiveNorm(const types::Matrix & stress, - const types::RVector & normal, - const types::RVector & tangent) { - AKANTU_DEBUG_IN(); - - Real normal_contrib, tangent_contrib; - - normal_stress.mul(stress, normal); - tangential_stress.mul(stress, tangent); - - normal_contrib = normal_stress.dot(normal); - tangent_contrib = tangential_stress.dot(tangent); - - if (normal_contrib < 0) normal_contrib = 0; - - AKANTU_DEBUG_OUT(); - - return std::sqrt(normal_contrib*normal_contrib - + tangent_contrib*tangent_contrib/beta/beta); -} - /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::computeTraction(const Vector & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Vector::iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Vector::iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); Vector::const_iterator normal_it = normal.begin(spatial_dimension); Vector::iteratortraction_end = tractions(el_type, ghost_type).end(spatial_dimension); Vector::iteratorsigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); Vector::iteratordelta_max_it = delta_max(el_type, ghost_type).begin(); Vector::iteratordelta_c_it = delta_c(el_type, ghost_type).begin(); Vector::iteratordamage_it = damage(el_type, ghost_type).begin(); /// compute scalars Real beta2_kappa2 = beta*beta/kappa/kappa; Real beta2_kappa = beta*beta/kappa; Real epsilon = std::numeric_limits::epsilon(); Real * memory_space = new Real[2*spatial_dimension]; types::Vector normal_opening(memory_space, spatial_dimension); types::Vector tangential_opening(memory_space + spatial_dimension, spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm; delta *= delta * beta2_kappa2; /// don't consider penetration contribution if (normal_opening_norm > 0) delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); /// full damage case or zero displacement case if (delta >= *delta_c_it || std::abs(delta) <= std::abs(delta) * epsilon) { /// set traction to zero (*traction_it).clear(); *damage_it = delta >= *delta_c_it; *delta_max_it = *damage_it * (*delta_c_it); } /// element not fully damaged else { /** * Compute traction @f$ \mathbf{T} = \left( * \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n * \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1- * \frac{\delta}{\delta_c} \right)@f$ */ *traction_it = tangential_opening; *traction_it *= beta2_kappa; /// update maximum displacement *delta_max_it = std::max(*delta_max_it, delta); *damage_it = *delta_max_it / *delta_c_it; Real k = *sigma_c_it / *delta_max_it * (1. - *damage_it); *traction_it *= k; /// use penalty coefficient in case of penetration if (normal_opening_norm < 0) k = penalty; normal_opening *= k; *traction_it += normal_opening; } } delete [] memory_space; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Material<_cohesive_linear> [" << std::endl; stream << space << " + sigma_c : " << sigma_c << std::endl; stream << space << " + beta : " << beta << std::endl; stream << space << " + G_cI : " << G_cI << std::endl; stream << space << " + G_cII : " << G_cII << std::endl; stream << space << " + rand : " << rand << std::endl; stream << space << " + penalty : " << penalty << std::endl; if(this->isInit()) { stream << space << " + kappa : " << kappa << std::endl; } MaterialCohesive::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialCohesiveLinearExtrinsic); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh index 761fc7cfd..93967aa8b 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh @@ -1,157 +1,157 @@ /** * @file material_cohesive_linear.hh * @author Marco Vocialta * @date Mon Feb 20 12:00:34 2012 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_EXTRINSIC_HH__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_EXTRINSIC_HH__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * Cohesive material linear damage for extrinsic case * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) * - beta : weighting parameter for sliding and normal opening (default: 0) * - G_cI : fracture energy for mode I (default: 0) * - G_cII : fracture energy for mode II (default: 0) * - rand : randomness factor (default: 0) * - penalty : stiffness in compression to prevent penetration */ template class MaterialCohesiveLinearExtrinsic : public MaterialCohesive { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialCohesiveLinearExtrinsic(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialCohesiveLinearExtrinsic(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set patameters virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// initialize the material computed parameter virtual void initMaterial(); /// resize vectors for new cohesive elements virtual void resizeCohesiveVectors(); /// compute stress norms on quadrature points for each facet for stress check virtual void computeStressNorms(const Vector & facet_stress, - types::RVector & stress_check); + Vector & stress_check); protected: /// constitutive law void computeTraction(const Vector & normal, ElementType el_type, GhostType ghost_type = _not_ghost); /// compute effective stress norm for insertion check inline Real computeEffectiveNorm(const types::Matrix & stress, const types::RVector & normal, const types::RVector & tangent); void computeTangentStiffness( __attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Vector & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } // void computeTangentStiffness(__attribute__((unused)) Vector & tangent_matrix, // __attribute__((unused)) const Vector & normal, // __attribute__((unused)) const ElementType & el_type, // __attribute__((unused)) GhostType ghost_type = _not_ghost) { // AKANTU_DEBUG_TO_IMPLEMENT(); // } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// critical effective stress ByElementTypeReal sigma_c_eff; /// beta parameter Real beta; /// mode I fracture energy Real G_cI; /// mode II fracture energy Real G_cII; /// kappa parameter Real kappa; /// penalty coefficient Real penalty; /// critical displacement ByElementTypeReal delta_c; /// vector to temporarily store the normal stress for the norm types::RVector normal_stress; /// vector to temporarily store the tangential stress for the norm types::RVector tangential_stress; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_cohesive_linear_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_EXTRINSIC_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc index b0752ace1..a49e5a548 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc @@ -1,633 +1,644 @@ /** * @file material_cohesive.cc * @author Marco Vocialta * @date Tue Feb 7 18:24:52 2012 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model,id), reversible_energy("reversible_energy", id), total_energy("total_energy", id), tractions_old("tractions (old)",id), opening_old("opening (old)",id), tractions("tractions",id), opening("opening",id), delta_max("delta max",id), damage("damage", id) { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); initInternalVector(reversible_energy, 1, _ek_cohesive); initInternalVector(total_energy, 1, _ek_cohesive); initInternalVector(tractions_old, spatial_dimension, _ek_cohesive); initInternalVector(tractions, spatial_dimension, _ek_cohesive); initInternalVector(opening_old, spatial_dimension, _ek_cohesive); initInternalVector(opening, spatial_dimension, _ek_cohesive); initInternalVector(delta_max, 1, _ek_cohesive); initInternalVector(damage, 1, _ek_cohesive); sigma_c = 0; rand = 0; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool MaterialCohesive::setParam(const std::string & key, const std::string & value, const ID & id) { return Material::setParam(key,value,id); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); fem_cohesive = &(model->getFEMClass("CohesiveFEM")); Mesh & mesh = fem_cohesive->getMesh(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, gt, _ek_cohesive); for(; it != last_type; ++it) { ElementType type = *it; element_filter.alloc(0, 1, type); } } resizeCohesiveVectors(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::resizeCohesiveVectors() { resizeInternalVector(reversible_energy, _ek_cohesive); resizeInternalVector(total_energy, _ek_cohesive); resizeInternalVector(tractions_old, _ek_cohesive); resizeInternalVector(tractions, _ek_cohesive); resizeInternalVector(opening_old, _ek_cohesive); resizeInternalVector(opening, _ek_cohesive); resizeInternalVector(delta_max, _ek_cohesive); resizeInternalVector(damage, _ek_cohesive); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::checkInsertion(const Vector & facet_stress, Vector & facet_insertion) { AKANTU_DEBUG_IN(); Vector & facets_check = model->getFacetsCheck(); ElementType type_facet = model->getFacetType(); UInt nb_quad_facet = model->getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_facet = facets_check.getSize(); - types::RVector stress_check(nb_quad_facet); + Vector stress_check(nb_facet, nb_quad_facet); + stress_check.clear(); computeStressNorms(facet_stress, stress_check); - for (UInt f = 0; f < nb_facet; ++f) { - if (facets_check(f) == true) { + bool * facet_check_it = facets_check.storage(); + Vector::iterator stress_check_it = + stress_check.begin(nb_quad_facet); + + Real * sigma_limit_it = model->getSigmaLimit().storage(); + + for (UInt f = 0; f < nb_facet; + ++f, ++facet_check_it, ++stress_check_it, ++sigma_limit_it) { + if (*facet_check_it == true) { for (UInt q = 0; q < nb_quad_facet; ++q) { - if (stress_check(q) > model->getSigmaLimit()(f)) { + if ((*stress_check_it)(q) > *sigma_limit_it) { facet_insertion.push_back(f); for (UInt qs = 0; qs < nb_quad_facet; ++qs) - sigma_insertion.push_back(stress_check(qs)); - facets_check(f) = false; + sigma_insertion.push_back((*stress_check_it)(qs)); + *facet_check_it = false; break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} t_e N_e dS @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::updateResidual(__attribute__((unused)) Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// compute traction computeTraction(ghost_type); /// update and assemble residual assembleResidual(ghost_type); /// compute energies computeEnergies(); /// update old values Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { tractions_old(*it, ghost_type).copy(tractions(*it, ghost_type)); opening_old(*it, ghost_type).copy(opening(*it, ghost_type)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Vector & residual = const_cast &>(model->getResidual()); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { const Vector & shapes = fem_cohesive->getShapes(*it, ghost_type); Vector & elem_filter = element_filter(*it, ghost_type); Vector & traction = tractions(*it, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(*it, ghost_type); UInt nb_element = elem_filter.getSize(); /// compute @f$t_i N_a@f$ Real * shapes_val = shapes.storage(); UInt * elem_filter_val = elem_filter.storage(); Vector * shapes_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_filtered_val = shapes_filtered->values; for (UInt el = 0; el < nb_element; ++el) { shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } shapes_filtered_val = shapes_filtered->values; // multiply traction by shapes - Vector traction_cpy(traction); - traction_cpy.extendComponentsInterlaced(size_of_shapes, spatial_dimension); + Vector * traction_cpy = new Vector(traction); + traction_cpy->extendComponentsInterlaced(size_of_shapes, spatial_dimension); - Real * traction_cpy_val = traction_cpy.storage(); + Real * traction_cpy_val = traction_cpy->storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { for (UInt n = 0; n < size_of_shapes; ++n,++shapes_filtered_val) { for (UInt i = 0; i < spatial_dimension; ++i) { *traction_cpy_val++ *= *shapes_filtered_val; } } } } delete shapes_filtered; /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ - Vector int_t_N(nb_element, spatial_dimension*size_of_shapes, - "int_t_N"); + Vector * int_t_N = new Vector(nb_element, spatial_dimension*size_of_shapes, + "int_t_N"); - fem_cohesive->integrate(traction_cpy, int_t_N, + fem_cohesive->integrate(*traction_cpy, *int_t_N, spatial_dimension*size_of_shapes, *it, ghost_type, &elem_filter); - int_t_N.extendComponentsInterlaced(2, int_t_N.getNbComponent() ); + delete traction_cpy; + + int_t_N->extendComponentsInterlaced(2, int_t_N->getNbComponent() ); - Real * int_t_N_val = int_t_N.storage(); + Real * int_t_N_val = int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < size_of_shapes*spatial_dimension; ++n) int_t_N_val[n] *= -1.; int_t_N_val += nb_nodes_per_element*spatial_dimension; } /// assemble - model->getFEMBoundary().assembleVector(int_t_N, residual, + model->getFEMBoundary().assembleVector(*int_t_N, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, &elem_filter, 1); + delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(Vector & current_position, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(*it, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); const Vector & shapes = fem_cohesive->getShapes(*it, ghost_type); Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt size_of_shapes = shapes.getNbComponent(); Real * shapes_val = shapes.storage(); UInt * elem_filter_val = elem_filter.storage(); Vector * shapes_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_filtered_val = shapes_filtered->values; for (UInt el = 0; el < nb_element; ++el) { shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } shapes_filtered_val = shapes_filtered->values; /** * compute A matrix @f$ \mathbf{A} = \left[\begin{array}{c c c c c c c c c c c c} * 1 & 0 & 0 & 0& 0 & 0 & -1& 0 & 0 &0 &0 &0 \\ * 0 &1& 0&0 &0 &0 &0 & -1& 0& 0 & 0 &0 \\ * 0 &0& 1&0 &0 &0 &0 & 0& -1& 0 & 0 &0 \\ * 0 &0& 0&1 &0 &0 &0 & 0& 0& -1 & 0 &0 \\ * 0 &0& 0&0 &1 &0 &0 & 0& 0& 0 & -1 &0 \\ * 0 &0& 0&0 &0 &1 &0 & 0& 0& 0 & 0 &-1 * \end{array} \right]@f$ **/ UInt size_of_A = spatial_dimension*size_of_shapes*spatial_dimension*nb_nodes_per_element; Real * A = new Real[size_of_A]; memset(A, 0, size_of_A*sizeof(Real)); for ( UInt i = 0; i < spatial_dimension*size_of_shapes; ++i) { A[ spatial_dimension * nb_nodes_per_element * i + i ] = 1; A[ spatial_dimension * nb_nodes_per_element * i + i + spatial_dimension * size_of_shapes ] = -1; } /// compute traction computeTraction(ghost_type); /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} @f$ Vector * tangent_stiffness_matrix = new Vector(nb_element*nb_quadrature_points, spatial_dimension* spatial_dimension, "tangent_stiffness_matrix"); computeTangentStiffness(*it, *tangent_stiffness_matrix, ghost_type); UInt size_of_N = spatial_dimension*size_of_shapes*spatial_dimension; Real * N = new Real[size_of_N]; UInt size_of_N_A = spatial_dimension*nb_nodes_per_element*spatial_dimension; Real * N_A = new Real[size_of_N_A]; Real * D_N_A = new Real[size_of_N_A]; UInt offset_At_Nt_D_N_A = spatial_dimension*nb_nodes_per_element* spatial_dimension*nb_nodes_per_element; Vector * at_nt_d_n_a = new Vector (nb_element*nb_quadrature_points, offset_At_Nt_D_N_A, "A^t*N^t*D*N*A"); Real * At_Nt_D_N_A = at_nt_d_n_a->storage(); Real * D = tangent_stiffness_matrix->storage(); /** * compute the N matrix @f$ \mathbf{N} = \begin{array}{cccccc} * N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) * \end{array} @f$ **/ for (UInt e = 0; e < nb_element; ++e) { Real * shapes_val = shapes_filtered_val + e * nb_quadrature_points * size_of_shapes; for (UInt q = 0; q < nb_quadrature_points; ++q) { memset(N, 0, size_of_N * sizeof(Real)); for (UInt i = 0; i < spatial_dimension ; ++i) { Real * Nvoigt_tmp = N + i * (size_of_shapes * spatial_dimension) + i; Real * Nregular = shapes_val + q * (size_of_shapes) ; for (UInt n = 0; n < size_of_shapes; ++n) { *Nvoigt_tmp = *Nregular; Nvoigt_tmp += spatial_dimension; Nregular ++; } } /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} {\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ Math::matrix_matrix (spatial_dimension, spatial_dimension*nb_nodes_per_element, size_of_shapes*spatial_dimension, N, A, N_A); Math::matrix_matrix (spatial_dimension, spatial_dimension*nb_nodes_per_element, spatial_dimension, D, N_A, D_N_A); Math::matrixt_matrix(spatial_dimension*nb_nodes_per_element, spatial_dimension*nb_nodes_per_element, spatial_dimension, N_A, D_N_A, At_Nt_D_N_A); At_Nt_D_N_A += offset_At_Nt_D_N_A; D += spatial_dimension * spatial_dimension; } } delete [] N; delete [] N_A; delete [] D_N_A; delete tangent_stiffness_matrix; Vector * K_e = new Vector(nb_element, offset_At_Nt_D_N_A , "K_e"); fem_cohesive->integrate(*at_nt_d_n_a, *K_e, offset_At_Nt_D_N_A, *it, ghost_type, &elem_filter); delete at_nt_d_n_a; model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, *it, ghost_type, &elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeTangentStiffness(const ElementType & el_type, Vector & tangent_matrix, GhostType ghost_type) { UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(el_type, ghost_type); Vector normal(nb_quadrature_points, spatial_dimension, "normal"); computeNormal(model->getCurrentPosition(), normal, el_type, ghost_type); computeTangentStiffness(el_type, tangent_matrix, normal, ghost_type); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quadrature_points = nb_element*fem_cohesive->getNbQuadraturePoints(*it, ghost_type); Vector normal(nb_quadrature_points, spatial_dimension, "normal"); /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(*it, ghost_type), *it, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normal, *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Vector & position, Vector & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); #define COMPUTE_NORMAL(type) \ fem_cohesive->getShapeFunctions(). \ computeNormalsOnControlPoints(position, \ normal, \ ghost_type, \ &(element_filter(type, ghost_type))) AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Vector & displacement, Vector & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); #define COMPUTE_OPENING(type) \ fem_cohesive->getShapeFunctions(). \ interpolateOnControlPoints(displacement, \ opening, \ spatial_dimension, \ ghost_type, \ &(element_filter(type, ghost_type))) AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeEnergies() { AKANTU_DEBUG_IN(); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); Real * memory_space = new Real[2*spatial_dimension]; types::RVector b(memory_space, spatial_dimension); types::RVector h(memory_space + spatial_dimension, spatial_dimension); for(; it != last_type; ++it) { Vector::iterator erev = reversible_energy(*it, _not_ghost).begin(); Vector::iterator etot = total_energy(*it, _not_ghost).begin(); Vector::iterator traction_it = tractions(*it, _not_ghost).begin(spatial_dimension); Vector::iterator traction_old_it = tractions_old(*it, _not_ghost).begin(spatial_dimension); Vector::iterator opening_it = opening(*it, _not_ghost).begin(spatial_dimension); Vector::iterator opening_old_it = opening_old(*it, _not_ghost).begin(spatial_dimension); Vector::iteratortraction_end = tractions(*it, _not_ghost).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration b = *opening_it; b -= *opening_old_it; h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } } delete [] memory_space; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate the dissipated energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { erev += fem_cohesive->integrate(reversible_energy(*it, _not_ghost), *it, _not_ghost, &element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate the dissipated energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Vector dissipated_energy(total_energy(*it, _not_ghost)); dissipated_energy -= reversible_energy(*it, _not_ghost); edis += fem_cohesive->integrate(dissipated_energy, *it, _not_ghost, &element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if (type == "reversible") return getReversibleEnergy(); else if (type == "dissipated") return getDissipatedEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void MaterialCohesive::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Material Cohesive [" << std::endl; Material::printself(stream, indent + 1); stream << space << "]" << std::endl; } __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh index 89b1c9495..6075dcddb 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh @@ -1,267 +1,267 @@ /** * @file material_cohesive.hh * @author Marco Vocialta * @date Tue Feb 7 17:50:23 2012 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" #include "fem_template.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_HH__ #define __AKANTU_MATERIAL_COHESIVE_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } __BEGIN_AKANTU__ class MaterialCohesive : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate< IntegratorCohesive, ShapeCohesive > MyFEMCohesiveType; public: MaterialCohesive(SolidMechanicsModel& model, const ID & id = ""); virtual ~MaterialCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// initialize the material computed parameter virtual void initMaterial(); /// resize vectors for new cohesive elements virtual void resizeCohesiveVectors(); /// compute the residual for this material virtual void updateResidual(Vector & current_position, GhostType ghost_type = _not_ghost); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(__attribute__((unused)) Real h, __attribute__((unused)) const Element & element = ElementNull) { AKANTU_DEBUG_TO_IMPLEMENT(); } // /// add an element to the local mesh filter // __aka_inline__ void addElement(const ElementType & type, // UInt element, // const GhostType & ghost_type); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// check stress for cohesive elements' insertion virtual void checkInsertion(const Vector & facet_stress, Vector & facet_insertion); /// interpolate stress on given positions for each element virtual void interpolateStress(__attribute__((unused)) const ElementType type, __attribute__((unused)) const Vector & coordinates, __attribute__((unused)) Vector & result) {}; protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } // virtual void computeTangentStiffness( Vector & tangent_matrix, // const Vector & normal, // const ElementType & el_type, // GhostType ghost_type = _not_ghost); // void computeTangentStiffness(const ElementType & el_type, // Vector & tangent_matrix, // GhostType ghost_type = _not_ghost){ // AKANTU_DEBUG_TO_IMPLEMENT(); // } void computeNormal(const Vector & position, Vector & normal, ElementType type, GhostType ghost_type); void computeOpening(const Vector & displacement, Vector & normal, ElementType type, GhostType ghost_type); template void computeNormal(const Vector & position, Vector & normal, GhostType ghost_type); /// assemble residual void assembleResidual(GhostType ghost_type = _not_ghost); /// assemble stiffness void assembleStiffnessMatrix(Vector & current_position, GhostType ghost_type); // template // void assembleStiffnessMatrix(Vector & current_position, // const ElementType & type, // GhostType ghost_type); /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// constitutive law virtual void computeTraction(const Vector & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// compute reversible and total energies by element void computeEnergies(); /// compute stress norms on quadrature points for each facet for stress check virtual void computeStressNorms(__attribute__((unused)) const Vector & facet_stress, - __attribute__((unused)) types::RVector & stress_check) { + __attribute__((unused)) Vector & stress_check) { AKANTU_DEBUG_TO_IMPLEMENT(); }; protected: /// compute the tangent stiffness matrix for an element type virtual void computeTangentStiffness(const ElementType & el_type, Vector & tangent_matrix, GhostType ghost_type = _not_ghost); virtual void computeTangentStiffness(const ElementType & el_type, Vector & tangent_matrix, const Vector & normal, GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening, Real); /// get the traction AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Traction, tractions, Real); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// get energy virtual Real getEnergy(std::string type); /// get sigma_c AKANTU_GET_MACRO(SigmaC, sigma_c, Real); /// get rand AKANTU_GET_MACRO(RandFactor, rand, Real); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// reversible energy by quadrature point ByElementTypeReal reversible_energy; /// total energy by quadrature point ByElementTypeReal total_energy; /// traction in all elements and quadrature points (previous time step) ByElementTypeReal tractions_old; /// opening in all elements and quadrature points (previous time step) ByElementTypeReal opening_old; protected: /// traction in all elements and quadrature points ByElementTypeReal tractions; /// opening in all elements and quadrature points ByElementTypeReal opening; /// Link to the cohesive fem object in the model MyFEMCohesiveType * fem_cohesive; /// critical stress Real sigma_c; /// randomness factor Real rand; /// vector to store stresses on facets for element insertions Vector sigma_insertion; /// maximum displacement ByElementTypeReal delta_max; /// damage ByElementTypeReal damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc index b6606557f..65a89b50f 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc @@ -1,1017 +1,1029 @@ /** * @file solid_mechanics_model_cohesive.cc * @author Marco Vocialta * @date Thu Apr 19 10:42:11 2012 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), mesh_facets(mesh.getSpatialDimension(), mesh.getNodes().getID(), id, memory_id), - elements_quad_facets("elements_quad_facets", id), stress_on_facet("stress_on_facet", id), facet_stress(0, spatial_dimension * spatial_dimension, "facet_stress"), facets_to_cohesive_el(0, 2, "facets_to_cohesive_el"), fragment_to_element("fragment_to_element", id), fragment_velocity(0, spatial_dimension, "fragment_velocity"), fragment_center(0, spatial_dimension, "fragment_center") { AKANTU_DEBUG_IN(); registerFEMObject("CohesiveFEM", mesh, spatial_dimension); MeshUtils::buildAllFacets(mesh, mesh_facets); /// assign cohesive type Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); for (; it != last; ++it) { const Vector & connectivity = mesh_facets.getConnectivity(*it); if (connectivity.getSize() != 0) { type_facet = *it; break; } } if (type_facet == _segment_2) type_cohesive = _cohesive_2d_4; else if (type_facet == _segment_3) type_cohesive = _cohesive_2d_6; mesh.addConnectivityType(type_cohesive); nb_fragment = 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFull(std::string material_file, AnalysisMethod method) { SolidMechanicsModel::initFull(material_file, method); if(material_file != "") { initCohesiveMaterial(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initCohesiveMaterial() { AKANTU_DEBUG_IN(); /// find the cohesive index cohesive_index = 0; while ((dynamic_cast(materials[cohesive_index]) == NULL) && cohesive_index <= materials.size()) ++cohesive_index; AKANTU_DEBUG_ASSERT(cohesive_index != materials.size(), "No cohesive materials in the material input file"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModelCohesive::initModel() { SolidMechanicsModel::initModel(); getFEM("CohesiveFEM").initShapeFunctions(_not_ghost); getFEM("CohesiveFEM").initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initExtrinsic() { AKANTU_DEBUG_IN(); /// assign sigma_c to each facet MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); UInt nb_facet = mesh_facets.getNbElement(type_facet); sigma_lim.resize(nb_facet); const Real sigma_c = mat_cohesive->getSigmaC(); const Real rand = mat_cohesive->getRandFactor(); std::srand(time(NULL)); if (facets_check.getSize() < 1) { const Vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet); facets_check.resize(nb_facet); for (UInt f = 0; f < nb_facet; ++f) { if (element_to_facet(f)(1) != ElementNull) { facets_check(f) = true; sigma_lim(f) = sigma_c * (1 + std::rand()/(Real)RAND_MAX * rand); } else { facets_check(f) = false; sigma_lim(f) = std::numeric_limits::max(); } } } /// compute normals on facets registerFEMObject("FacetsFEM", mesh_facets, spatial_dimension-1); getFEM("FacetsFEM").initShapeFunctions(); getFEM("FacetsFEM").computeNormalsOnControlPoints(); /// THIS HAS TO BE CHANGED: /* ------------------------------------------------------------------------ */ const Vector & normals = getFEM("FacetsFEM").getNormalsOnQuadPoints(type_facet); Vector::const_iterator normal_it = normals.begin(spatial_dimension); tangents.resize(normals.getSize()); tangents.extendComponentsInterlaced(spatial_dimension, tangents.getNbComponent()); Vector::iterator tangent_it = tangents.begin(spatial_dimension); for (UInt i = 0; i < normals.getSize(); ++i, ++normal_it, ++tangent_it) { Math::normal2( (*normal_it).storage(), (*tangent_it).storage() ); } /* ------------------------------------------------------------------------ */ /// compute quadrature points coordinates on facets Vector & position = mesh.getNodes(); UInt nb_quad_per_facet = getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_tot_quad = nb_quad_per_facet * nb_facet; Vector quad_facets(nb_tot_quad, spatial_dimension); getFEM("FacetsFEM").interpolateOnQuadraturePoints(position, quad_facets, spatial_dimension, type_facet); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure + ByElementTypeReal elements_quad_facets; + Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element == 0) continue; /// compute elements' quadrature points and list of facet /// quadrature points positions by element Vector & facet_to_element = mesh_facets.getSubelementToElement(*it); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); elements_quad_facets.alloc(nb_element * nb_facet_per_elem * nb_quad_per_facet, spatial_dimension, *it); Vector & el_q_facet = elements_quad_facets(*it); stress_on_facet.alloc(el_q_facet.getSize(), spatial_dimension * spatial_dimension, *it); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { UInt global_facet = facet_to_element(el, f).element; for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < spatial_dimension; ++s) { el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s) = quad_facets(global_facet * nb_quad_per_facet + q, s); } } } } } /// initialize the interpolation function - materials[0]->initInterpolateElementalField(); + materials[0]->initElementalFieldInterpolation(elements_quad_facets); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); Vector facet_insertion; UInt nb_materials = materials.size(); UInt nb_quad_per_facet = getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_facet = mesh_facets.getNbElement(type_facet); /// vector containing stresses coming from the two elements of each facets facet_stress.resize(2 * nb_facet * nb_quad_per_facet); facet_stress_count.resize(nb_facet); facet_stress_count.clear(); /// loop over materials for (UInt m = 0; m < nb_materials; ++m) { if (m == cohesive_index) continue; Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); /// loop over element type for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element == 0) continue; - const Vector & el_q_facet = elements_quad_facets(*it); Vector & stress_on_f = stress_on_facet(*it); /// interpolate stress on facet quadrature points positions materials[m]->interpolateStress(*it, - el_q_facet, stress_on_f); /// store the interpolated stresses on the facet_stress vector Vector & facet_to_element = mesh_facets.getSubelementToElement(*it); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); - for (UInt el = 0; el < nb_element; ++el) { + Vector::iterator > facet_to_el_it = + facet_to_element.begin(nb_facet_per_elem); + + Vector::iterator stress_on_f_it = + stress_on_f.begin(spatial_dimension, spatial_dimension); + + UInt sp2 = spatial_dimension * spatial_dimension; + UInt nb_quad_f_two = nb_quad_per_facet * 2; + + for (UInt el = 0; el < nb_element; ++el, ++facet_to_el_it) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { - UInt global_facet = facet_to_element(el, f).element; - - for (UInt q = 0; q < nb_quad_per_facet; ++q) { - for (UInt s = 0; s < spatial_dimension * spatial_dimension; ++s) { - facet_stress(2 * global_facet * nb_quad_per_facet - + facet_stress_count(global_facet) * nb_quad_per_facet - + q, s) - = stress_on_f(el * nb_facet_per_elem * nb_quad_per_facet - + f * nb_quad_per_facet + q, s); - } + UInt global_facet = (*facet_to_el_it)(f).element; + + for (UInt q = 0; q < nb_quad_per_facet; ++q, ++stress_on_f_it) { + types::Matrix facet_stress_local(facet_stress.storage() + + (global_facet * nb_quad_f_two + + q * 2 + + facet_stress_count(global_facet)) + * sp2, + spatial_dimension, + spatial_dimension); + + facet_stress_local = *stress_on_f_it; } - facet_stress_count(global_facet) = true; } } } } MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); mat_cohesive->checkInsertion(facet_stress, facet_insertion); if (facet_insertion.getSize() != 0) insertCohesiveElements(facet_insertion); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertCohesiveElements(const Vector & facet_insertion) { AKANTU_DEBUG_IN(); if(facet_insertion.getSize() == 0) return; Vector facet_material(facet_insertion.getSize(), 1, cohesive_index); insertCohesiveElements(facet_insertion, facet_material); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertCohesiveElements(const Vector & facet_insertion, const Vector & facet_material) { AKANTU_DEBUG_IN(); if(facet_insertion.getSize() == 0) return; Vector doubled_nodes(0, 2); Vector doubled_facets(0, 2); /// update mesh for (UInt f = 0; f < facet_insertion.getSize(); ++f) { Element facet(type_facet, facet_insertion(f)); doubleFacet(facet, doubled_nodes, doubled_facets); } /// double middle nodes if it's the case if (type_facet == _segment_3) doubleMiddleNode(doubled_nodes, doubled_facets); /// loop over doubled facets to insert cohesive elements Vector & conn_cohesive = mesh.getConnectivity(type_cohesive); const Vector & conn_facet = mesh_facets.getConnectivity(type_facet); UInt nb_nodes_per_facet = conn_facet.getNbComponent(); const Vector & position = mesh.getNodes(); Vector & element_mat = getElementMaterial(type_cohesive); Vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet); const Real epsilon = std::numeric_limits::epsilon(); for (UInt f = 0; f < doubled_facets.getSize(); ++f) { UInt nb_cohesive_elements = conn_cohesive.getSize(); conn_cohesive.resize(nb_cohesive_elements + 1); UInt first_facet = doubled_facets(f, 0); UInt second_facet = doubled_facets(f, 1); /// copy first facet's connectivity for (UInt n = 0; n < nb_nodes_per_facet; ++n) conn_cohesive(nb_cohesive_elements, n) = conn_facet(first_facet, n); /// check if first nodes of the two facets are coincident or not UInt first_facet_node = conn_facet(first_facet, 0); UInt second_facet_node = conn_facet(second_facet, 0); bool second_facet_inversion = false; for (UInt dim = 0; dim < mesh.getSpatialDimension(); ++dim) { if (std::abs( (position(first_facet_node, dim) - position(second_facet_node, dim)) / position(second_facet_node, dim)) >= epsilon) { second_facet_inversion = true; break; } } /// if the two nodes are coincident second facet connectivity is /// normally copied, otherwise the copy is reverted if (!second_facet_inversion) { for (UInt n = 0; n < nb_nodes_per_facet; ++n) conn_cohesive(nb_cohesive_elements, n + nb_nodes_per_facet) = conn_facet(second_facet, n); } else { for (UInt n = 0; n < nb_nodes_per_facet; ++n) conn_cohesive(nb_cohesive_elements, n + nb_nodes_per_facet) = conn_facet(second_facet, nb_nodes_per_facet - n - 1); } /// assign the cohesive material to the new element element_mat.resize(nb_cohesive_elements + 1); element_mat(nb_cohesive_elements) = facet_material(f); materials[cohesive_index]->addElement(type_cohesive, nb_cohesive_elements, _not_ghost); /// update element_to_facet vectors Element cohesive_element(type_cohesive, nb_cohesive_elements, _not_ghost, _ek_cohesive); element_to_facet(first_facet)(1) = cohesive_element; element_to_facet(second_facet)(1) = cohesive_element; /// update facets_to_cohesive_el vector facets_to_cohesive_el.resize(nb_cohesive_elements + 1); facets_to_cohesive_el(nb_cohesive_elements, 0) = first_facet; facets_to_cohesive_el(nb_cohesive_elements, 1) = second_facet; } /// update shape functions getFEM("CohesiveFEM").initShapeFunctions(_not_ghost); /// resize cohesive material vectors MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); AKANTU_DEBUG_ASSERT(mat_cohesive, "No cohesive materials in the materials vector"); mat_cohesive->resizeCohesiveVectors(); /// update nodal values updateDoubledNodes(doubled_nodes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::doubleMiddleNode(Vector & doubled_nodes, const Vector & doubled_facets) { AKANTU_DEBUG_IN(); Vector & conn_facet = mesh_facets.getConnectivity(type_facet); Vector & position = mesh.getNodes(); Vector > & elem_to_facet = mesh_facets.getElementToSubelement(type_facet); for (UInt f = 0; f < doubled_facets.getSize(); ++f) { UInt facet_first = doubled_facets(f, 0); UInt facet_second = doubled_facets(f, 1); UInt new_node = position.getSize(); /// store doubled nodes UInt nb_doubled_nodes = doubled_nodes.getSize(); doubled_nodes.resize(nb_doubled_nodes + 1); UInt old_node = conn_facet(facet_first, 2); doubled_nodes(nb_doubled_nodes, 0) = old_node; doubled_nodes(nb_doubled_nodes, 1) = new_node; /// update position position.resize(new_node + 1); for (UInt dim = 0; dim < spatial_dimension; ++dim) position(new_node, dim) = position(old_node, dim); /// update facet connectivity conn_facet(facet_second, 2) = new_node; /// update element connectivity for (UInt el = 0; el < elem_to_facet(facet_second).getSize(); ++el) { const ElementType type_elem = elem_to_facet(facet_second)(el).type; if (type_elem != _not_defined) { UInt elem_global = elem_to_facet(facet_second)(el).element; Vector & conn_elem = mesh.getConnectivity(type_elem); for (UInt n = 0; n < conn_elem.getNbComponent(); ++n) { if (conn_elem(elem_global, n) == old_node) conn_elem(elem_global, n) = new_node; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateDoubledNodes(const Vector & doubled_nodes) { AKANTU_DEBUG_IN(); UInt nb_old_nodes = displacement->getSize(); UInt nb_new_nodes = nb_old_nodes + doubled_nodes.getSize(); displacement ->resize(nb_new_nodes); velocity ->resize(nb_new_nodes); acceleration ->resize(nb_new_nodes); if (increment_acceleration) increment_acceleration->resize(nb_new_nodes); force ->resize(nb_new_nodes); residual ->resize(nb_new_nodes); boundary ->resize(nb_new_nodes); mass ->resize(nb_new_nodes); if (increment) increment->resize(nb_new_nodes); //initsolver(); /** * @todo temporary patch, should be done in a better way that works * always (pbc, parallel, ...) **/ delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); if (method != _explicit_dynamic) { if(method == _implicit_dynamic) delete stiffness_matrix; delete jacobian_matrix; delete solver; SolverOptions solver_options; initImplicit((method == _implicit_dynamic), solver_options); } for (UInt n = 0; n < doubled_nodes.getSize(); ++n) { UInt old_node = doubled_nodes(n, 0); UInt new_node = doubled_nodes(n, 1); for (UInt dim = 0; dim < displacement->getNbComponent(); ++dim) { (*displacement)(new_node, dim) = (*displacement)(old_node, dim); } for (UInt dim = 0; dim < velocity->getNbComponent(); ++dim) { (*velocity)(new_node, dim) = (*velocity)(old_node, dim); } for (UInt dim = 0; dim < acceleration->getNbComponent(); ++dim) { (*acceleration)(new_node, dim) = (*acceleration)(old_node, dim); } if (increment_acceleration) { for (UInt dim = 0; dim < increment_acceleration->getNbComponent(); ++dim) { (*increment_acceleration)(new_node, dim) = (*increment_acceleration)(old_node, dim); } } } if (increment) { for (UInt n = 0; n < doubled_nodes.getSize(); ++n) { UInt old_node = doubled_nodes(n, 0); UInt new_node = doubled_nodes(n, 1); for (UInt dim = 0; dim < increment->getNbComponent(); ++dim) (*increment)(new_node, dim) = (*increment)(old_node, dim); } } assembleMassLumped(); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::doubleFacet(Element & facet, Vector & doubled_nodes, Vector & doubled_facets) { AKANTU_DEBUG_IN(); const UInt f_index = facet.element; const ElementType type_facet = facet.type; const ElementType type_subfacet = mesh.getFacetElementType(type_facet); const UInt nb_subfacet = mesh.getNbFacetsPerElement(type_facet); Vector > & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet); Vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet); Vector & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet); /// adding a new facet by copying original one /// create new connectivity Vector & conn_facet = mesh_facets.getConnectivity(type_facet); UInt nb_facet = conn_facet.getSize(); conn_facet.resize(nb_facet + 1); for (UInt n = 0; n < conn_facet.getNbComponent(); ++n) conn_facet(nb_facet, n) = conn_facet(f_index, n); /// store doubled facets UInt nb_doubled_facets = doubled_facets.getSize(); doubled_facets.resize(nb_doubled_facets + 1); doubled_facets(nb_doubled_facets, 0) = f_index; doubled_facets(nb_doubled_facets, 1) = nb_facet; /// update elements connected to facet Vector first_facet_list = element_to_facet(f_index); element_to_facet.push_back(first_facet_list); /// set new and original facets as boundary facets + AKANTU_DEBUG_ASSERT(element_to_facet(f_index)(1) != ElementNull, + "can't double a facet on the boundary!"); + element_to_facet(f_index)(1) = ElementNull; element_to_facet(nb_facet)(0) = element_to_facet(nb_facet)(1); element_to_facet(nb_facet)(1) = ElementNull; /// update facet_to_element vector ElementType type = element_to_facet(nb_facet)(0).type; UInt el = element_to_facet(nb_facet)(0).element; Vector & facet_to_element = mesh_facets.getSubelementToElement(type); UInt i; for (i = 0; facet_to_element(el, i).element != f_index && i <= facet_to_element.getNbComponent(); ++i); facet_to_element(el, i).element = nb_facet; /// create new links to subfacets and update list of facets /// connected to subfacets subfacet_to_facet.resize(nb_facet + 1); for (UInt sf = 0; sf < nb_subfacet; ++sf) { subfacet_to_facet(nb_facet, sf) = subfacet_to_facet(f_index, sf); UInt sf_index = subfacet_to_facet(f_index, sf).element; /// find index to start looping around facets connected to current /// subfacet UInt start = 0; UInt nb_connected_facets = facet_to_subfacet(sf_index).getSize(); while (facet_to_subfacet(sf_index)(start).element != f_index && start <= facet_to_subfacet(sf_index).getSize()) ++start; /// add the new facet to the list next to the original one ++nb_connected_facets; facet_to_subfacet(sf_index).resize(nb_connected_facets); for (UInt f = nb_connected_facets - 1; f > start; --f) { facet_to_subfacet(sf_index)(f) = facet_to_subfacet(sf_index)(f - 1); } /// check if the new facet should be inserted before or after the /// original one: the second element connected to both original /// and new facet will be _not_defined, so I check if the first /// one is equal to one of the elements connected to the following /// facet in the facet_to_subfacet vector UInt f_start = facet_to_subfacet(sf_index)(start).element; UInt f_next; if (start + 2 == nb_connected_facets) f_next = facet_to_subfacet(sf_index)(0).element; else f_next = facet_to_subfacet(sf_index)(start + 2).element; if ((element_to_facet(f_start)(0) == element_to_facet(f_next)(0)) || ( element_to_facet(f_start)(0) == element_to_facet(f_next)(1))) facet_to_subfacet(sf_index)(start).element = nb_facet; else facet_to_subfacet(sf_index)(start + 1).element = nb_facet; /// loop on every facet connected to the current subfacet for (UInt f = start + 2; ; ++f) { /// reset f in order to continue looping from the beginning if (f == nb_connected_facets) f = 0; /// exit loop if it reaches the end if (f == start) break; /// if current loop facet is on the boundary, double subfacet UInt f_global = facet_to_subfacet(sf_index)(f).element; if (element_to_facet(f_global)(1).type == _not_defined || element_to_facet(f_global)(1).kind == _ek_cohesive) { doubleSubfacet(subfacet_to_facet(f_index, sf), start, f, doubled_nodes); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::doubleSubfacet(const Element & subfacet, UInt start, UInt end, Vector & doubled_nodes) { AKANTU_DEBUG_IN(); const UInt sf_index = subfacet.element; const ElementType type_subfacet = subfacet.type; Vector > & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet); UInt nb_subfacet = facet_to_subfacet.getSize(); Vector & conn_point = mesh_facets.getConnectivity(_point); Vector & position = mesh.getNodes(); /// add the new subfacet if (spatial_dimension == 2) { UInt new_node = position.getSize(); /// add new node in connectivity UInt new_subfacet = conn_point.getSize(); conn_point.resize(new_subfacet + 1); conn_point(new_subfacet) = new_node; /// store doubled nodes UInt nb_doubled_nodes = doubled_nodes.getSize(); doubled_nodes.resize(nb_doubled_nodes + 1); UInt old_node = doubled_nodes(nb_doubled_nodes, 0) = conn_point(sf_index); doubled_nodes(nb_doubled_nodes, 1) = new_node; /// update position position.resize(new_node + 1); for (UInt dim = 0; dim < spatial_dimension; ++dim) position(new_node, dim) = position(old_node, dim); } /// create a vector for the new subfacet in facet_to_subfacet facet_to_subfacet.resize(nb_subfacet + 1); UInt nb_connected_facets = facet_to_subfacet(sf_index).getSize(); /// loop over facets from start to end for (UInt f = start + 1; ; ++f) { /// reset f in order to continue looping from the beginning if (f == nb_connected_facets) f = 0; UInt f_global = facet_to_subfacet(sf_index)(f).element; ElementType type_facet = facet_to_subfacet(sf_index)(f).type; Vector & conn_facet = mesh_facets.getConnectivity(type_facet); UInt old_node = conn_point(sf_index); UInt new_node = conn_point(conn_point.getSize() - 1); /// update facet connectivity UInt i; for (i = 0; conn_facet(f_global, i) != old_node && i <= conn_facet.getNbComponent(); ++i); conn_facet(f_global, i) = new_node; /// update element connectivity Vector > & elem_to_facet = mesh_facets.getElementToSubelement(type_facet); for (UInt el = 0; el < elem_to_facet(f_global).getSize(); ++el) { const ElementType type_elem = elem_to_facet(f_global)(el).type; if (type_elem != _not_defined) { UInt elem_global = elem_to_facet(f_global)(el).element; Vector & conn_elem = mesh.getConnectivity(type_elem); for (UInt n = 0; n < conn_elem.getNbComponent(); ++n) { if (conn_elem(elem_global, n) == old_node) conn_elem(elem_global, n) = new_node; } } } /// update subfacet_to_facet vector Vector & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet); for (i = 0; subfacet_to_facet(f_global, i).element != sf_index && i <= subfacet_to_facet.getNbComponent(); ++i); subfacet_to_facet(f_global, i).element = nb_subfacet; /// add current facet to facet_to_subfacet last position Element current_facet(type_facet, f_global); facet_to_subfacet(nb_subfacet).push_back(current_facet); /// exit loop if it reaches the end if (f == end) break; } /// rearrange the facets connected to the original subfacet and /// compute the new number of facets connected to it if (end < start) { for (UInt f = 0; f < start - end; ++f) facet_to_subfacet(sf_index)(f) = facet_to_subfacet(sf_index)(f + end + 1); nb_connected_facets = start - end; } else { for (UInt f = 1; f < nb_connected_facets - end; ++f) facet_to_subfacet(sf_index)(start + f) = facet_to_subfacet(sf_index)(end + f); nb_connected_facets -= end - start; } /// resize list of facets of the original subfacet facet_to_subfacet(sf_index).resize(nb_connected_facets); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::buildFragmentsList() { AKANTU_DEBUG_IN(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); const UInt max = std::numeric_limits::max(); for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element != 0) { /// initialize the list of checked elements and the list of /// elements to be checked fragment_to_element.alloc(nb_element, 1, *it); Vector & frag_to_el = fragment_to_element(*it); for (UInt el = 0; el < nb_element; ++el) frag_to_el(el) = max; } } Vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet); nb_fragment = 0; it = mesh.firstType(spatial_dimension); MaterialCohesive * mat_cohesive = dynamic_cast(materials[cohesive_index]); const Vector & damage = mat_cohesive->getDamage(type_cohesive); UInt nb_quad_cohesive = getFEM("CohesiveFEM").getNbQuadraturePoints(type_cohesive); Real epsilon = std::numeric_limits::epsilon(); for (; it != last; ++it) { Vector & checked_el = fragment_to_element(*it); UInt nb_element = checked_el.getSize(); if (nb_element != 0) { Vector & facet_to_element = mesh_facets.getSubelementToElement(*it); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); /// loop on elements for (UInt el = 0; el < nb_element; ++el) { if (checked_el(el) == max) { /// build fragment ++nb_fragment; checked_el(el) = nb_fragment - 1; Vector elem_to_check; Element first_el(*it, el); elem_to_check.push_back(first_el); /// keep looping while there are elements to check while (elem_to_check.getSize() != 0) { UInt nb_elem_check = elem_to_check.getSize(); for (UInt el_check = 0; el_check < nb_elem_check; ++el_check) { Element current_el = elem_to_check(el_check); for (UInt f = 0; f < nb_facet_per_elem; ++f) { /// find adjacent element on current facet UInt global_facet = facet_to_element(current_el.element, f).element; Element next_el; for (UInt i = 0; i < 2; ++i) { next_el = element_to_facet(global_facet)(i); if (next_el != current_el) break; } if (next_el.kind == _ek_cohesive) { /// fragmention occurs when the cohesive element has /// reached damage = 1 on every quadrature point UInt q = 0; while (q < nb_quad_cohesive && std::abs(damage(next_el.element * nb_quad_cohesive + q) - 1) <= epsilon) ++q; if (q == nb_quad_cohesive) next_el = ElementNull; else { /// check which facet is the correct one UInt other_facet_index = facets_to_cohesive_el(next_el.element, 0) == global_facet; UInt other_facet = facets_to_cohesive_el(next_el.element, other_facet_index); /// get the other regualar element next_el = element_to_facet(other_facet)(0); } } /// if it exists, add it to the fragment list if (next_el != ElementNull) { Vector & checked_next_el = fragment_to_element(next_el.type); /// check if the element isn't already part of a fragment if (checked_next_el(next_el.element) == max) { checked_next_el(next_el.element) = nb_fragment - 1; elem_to_check.push_back(next_el); } } } } /// erase elements that have already been checked for (UInt el_check = nb_elem_check; el_check > 0; --el_check) elem_to_check.erase(el_check - 1); } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeFragmentsMV() { AKANTU_DEBUG_IN(); fragment_mass.resize(nb_fragment); fragment_mass.clear(); fragment_velocity.resize(nb_fragment); fragment_velocity.clear(); fragment_center.resize(nb_fragment); fragment_center.clear(); UInt nb_nodes = mesh.getNbNodes(); Vector checked_node(nb_nodes); checked_node.clear(); const Vector & position = mesh.getNodes(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element != 0) { const Vector & frag_to_el = fragment_to_element(*it); const Vector & connectivity = mesh.getConnectivity(*it); UInt nb_nodes_per_elem = connectivity.getNbComponent(); /// loop over each node of every element for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = connectivity(el, n); /// if the node hasn't been checked, store its data if (!checked_node(node)) { fragment_mass(frag_to_el(el)) += (*mass)(node); for (UInt s = 0; s < spatial_dimension; ++s) { fragment_velocity(frag_to_el(el), s) += (*mass)(node) * (*velocity)(node, s); fragment_center(frag_to_el(el), s) += (*mass)(node) * position(node, s); } checked_node(node) = true; } } } } } for (UInt frag = 0; frag < nb_fragment; ++frag) { for (UInt s = 0; s < spatial_dimension; ++s) { fragment_velocity(frag, s) /= fragment_mass(frag); fragment_center(frag, s) /= fragment_mass(frag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeFragmentsData() { AKANTU_DEBUG_IN(); buildFragmentsList(); computeFragmentsMV(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "SolidMechanicsModelCohesive [" << std::endl; SolidMechanicsModel::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh index c9c60461c..3fb2e0832 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh @@ -1,237 +1,234 @@ /** * @file solid_mechanics_model_cohesive.hh * @author Marco Vocialta * @date Thu Apr 19 10:34:00 2012 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ __BEGIN_AKANTU__ class SolidMechanicsModelCohesive : public SolidMechanicsModel { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate< IntegratorCohesive, ShapeCohesive > MyFEMCohesiveType; SolidMechanicsModelCohesive(Mesh & mesh, UInt spatial_dimension = 0, const ID & id = "solid_mechanics_model_cohesive", const MemoryID & memory_id = 0); // virtual ~SolidMechanicsModelCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// function to perform a stress check on each facet and insert /// cohesive elements if needed void checkCohesiveStress(); /// function to insert cohesive elements on the selected facets void insertCohesiveElements(const Vector & facet_insertion); /// function to insert cohesive elements on the selected facets by specifying /// the material void insertCohesiveElements(const Vector & facet_insertion, const Vector & facet_material); /// initialize completely the model void initFull(std::string material_file, AnalysisMethod method = _explicit_dynamic); /// initialize completely the model for extrinsic elements void initExtrinsic(); /// initialize the model void initModel(); /// initialize cohesive material void initCohesiveMaterial(); /// build fragments and compute their data (mass, velocity..) void computeFragmentsData(); private: /// function to double a given facet and update the list of doubled /// nodes void doubleFacet(Element & facet, Vector & doubled_nodes, Vector & doubled_facets); /// function to double a subfacet given start and end index for /// local facet_to_subfacet vector, and update the list of doubled /// nodes void doubleSubfacet(const Element & subfacet, UInt start, UInt end, Vector & doubled_nodes); /// double middle nodes if facets are _segment_3 void doubleMiddleNode(Vector & doubled_nodes, const Vector & doubled_facets); /// function to update nodal parameters for doubled nodes void updateDoubledNodes(const Vector & doubled_nodes); /// build fragments list void buildFragmentsList(); /// compute fragments' mass and velocity void computeFragmentsMV(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the cohesive element type AKANTU_GET_MACRO(CohesiveElementType, type_cohesive, ElementType); /// get the cohesive index AKANTU_GET_MACRO(CohesiveIndex, cohesive_index, UInt); /// get the facet type AKANTU_GET_MACRO(FacetType, type_facet, ElementType); /// get the facet mesh AKANTU_GET_MACRO(MeshFacets, mesh_facets, const Mesh &); /// get the sigma limit vector for automatic insertion AKANTU_GET_MACRO(SigmaLimit, sigma_lim, const Vector &); AKANTU_GET_MACRO_NOT_CONST(SigmaLimit, sigma_lim, Vector &); /// get the facets check vector AKANTU_GET_MACRO_NOT_CONST(FacetsCheck, facets_check, Vector &); /// get the stress on facets vector AKANTU_GET_MACRO(StressOnFacets, facet_stress, const Vector &); /// get the number of fragments AKANTU_GET_MACRO(NbFragment, nb_fragment, UInt); /// get the fragment_to_element vectors AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FragmentToElement, fragment_to_element, UInt); /// get mass for each fragment AKANTU_GET_MACRO(FragmentsMass, fragment_mass, const Vector &); /// get average velocity for each fragment AKANTU_GET_MACRO(FragmentsVelocity, fragment_velocity, const Vector &); /// get the center of mass coordinates for each fragment AKANTU_GET_MACRO(FragmentsCenter, fragment_center, const Vector &); /// THIS HAS TO BE CHANGED AKANTU_GET_MACRO(Tangents, tangents, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// cohesive element type ElementType type_cohesive; /// facet type ElementType type_facet; /// cohesive material index in materials vector UInt cohesive_index; /// mesh containing facets and their data structures Mesh mesh_facets; /// vector containing a sigma limit for automatic insertion Vector sigma_lim; /// vector containing facets in which cohesive elements can be automatically inserted Vector facets_check; /// @todo store tangents when normals are computed: Vector tangents; - /// list of facet quadrature points positions by element - ByElementTypeReal elements_quad_facets; - /// list of stresses on facet quadrature points for every element ByElementTypeReal stress_on_facet; /// already counted facets in stress check Vector facet_stress_count; /// stress on facets on the two sides by quadrature point Vector facet_stress; /// list of facets connected to each cohesive element Vector facets_to_cohesive_el; /// fragment number for each element ByElementTypeUInt fragment_to_element; /// number of fragments UInt nb_fragment; /// mass for each fragment Vector fragment_mass; /// average velocity for each fragment Vector fragment_velocity; /// center of mass coordinates for each element Vector fragment_center; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "solid_mechanics_model_cohesive_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModelCohesive & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/mesh.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/mesh.geo index 4526cbf34..b7d030c25 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/mesh.geo +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/mesh.geo @@ -1,22 +1,13 @@ h = 1.0; -Point(1) = {0, 0, 0, h}; -Point(2) = {0, 1, 0, h}; -Point(3) = {-1, 1, 0, h}; -Point(4) = {-1, 0, 0, h}; -Point(6) = {1, 0, 0, h}; -Point(7) = {1, 1, 0, h}; -Point(8) = {1, -1, 0, h}; -Point(9) = {0, -1, 0, h}; -Point(10) = {-1, -1, 0, h}; -Line(1) = {10, 9}; -Line(2) = {9, 8}; -Line(3) = {8, 6}; -Line(4) = {6, 7}; -Line(5) = {7, 2}; -Line(6) = {2, 3}; -Line(7) = {3, 4}; -Line(8) = {4, 10}; -Line Loop(9) = {3, 4, 5, 6, 7, 8, 1, 2}; -Plane Surface(10) = {9}; -Physical Surface(11) = {10}; +Point(1) = {1, 1, 0, h}; +Point(2) = {-1, 1, 0, h}; +Point(3) = {-1, -1, 0, h}; +Point(4) = {1, -1, 0, h}; +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; +Line Loop(5) = {4, 1, 2, 3}; +Plane Surface(6) = {5}; +Physical Surface(7) = {6}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc index afc161753..de86ded56 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc @@ -1,250 +1,250 @@ /** * @file test_cohesive.cc * @author Marco Vocialta * @date Fri Feb 24 14:32:31 2012 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model_cohesive.hh" #include "material.hh" #include "io_helper.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char *argv[]) { initialize(argc, argv); debug::setDebugLevel(dblWarning); const UInt spatial_dimension = 2; const UInt max_steps = 1000; const ElementType type = _triangle_6; Mesh mesh(spatial_dimension); MeshIOMSH mesh_io; mesh_io.read("mesh.msh", mesh); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull("material.dat"); model.initExtrinsic(); Real time_step = model.getStableTimeStep()*0.05; model.setTimeStep(time_step); // std::cout << "Time step: " << time_step << std::endl; model.assembleMassLumped(); /* ------------------------------------------------------------------------ */ /* Facet part */ /* ------------------------------------------------------------------------ */ // std::cout << mesh << std::endl; const Mesh & mesh_facets = model.getMeshFacets(); const ElementType type_facet = mesh.getFacetElementType(type); UInt nb_facet = mesh_facets.getNbElement(type_facet); const Vector & position = mesh.getNodes(); // const Vector & connectivity = mesh_facets.getConnectivity(type_facet); Vector & sigma_lim = model.getSigmaLimit(); Vector & facet_check = model.getFacetsCheck(); Real * bary_facet = new Real[spatial_dimension]; for (UInt f = 0; f < nb_facet; ++f) { mesh_facets.getBarycenter(f, type_facet, bary_facet); if (bary_facet[1] < 0.30 && bary_facet[1] > 0.20) { sigma_lim(f) = 100; facet_check(f) = true; std::cout << f << std::endl; } else { sigma_lim(f) = 1e10; facet_check(f) = false; } } delete[] bary_facet; // std::cout << mesh << std::endl; /* ------------------------------------------------------------------------ */ /* End of facet part */ /* ------------------------------------------------------------------------ */ Vector & velocity = model.getVelocity(); Vector & boundary = model.getBoundary(); Vector & displacement = model.getDisplacement(); // const Vector & residual = model.getResidual(); UInt nb_nodes = mesh.getNbNodes(); /// boundary conditions for (UInt n = 0; n < nb_nodes; ++n) { if (position(n, 1) > 0.99 || position(n, 1) < -0.99) boundary(n, 1) = true; if (position(n, 0) > 0.99 || position(n, 0) < -0.99) boundary(n, 0) = true; } model.updateResidual(); - // iohelper::ElemType paraview_type = iohelper::TRIANGLE2; - // UInt nb_element = mesh.getNbElement(type); - - // /// initialize the paraview output - // iohelper::DumperParaview dumper; - // dumper.SetMode(iohelper::TEXT); - // dumper.SetPoints(mesh.getNodes().values, - // spatial_dimension, mesh.getNbNodes(), "explicit"); - // dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, - // paraview_type, nb_element, iohelper::C_MODE); - // dumper.AddNodeDataField(model.getDisplacement().values, - // spatial_dimension, "displacements"); - // dumper.AddNodeDataField(model.getVelocity().values, - // spatial_dimension, "velocity"); - // dumper.AddNodeDataField(model.getAcceleration().values, - // spatial_dimension, "acceleration"); - // dumper.AddNodeDataField(model.getForce().values, - // spatial_dimension, "applied_force"); - // dumper.AddNodeDataField(model.getResidual().values, - // spatial_dimension, "forces"); - // dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, - // spatial_dimension*spatial_dimension, "strain"); - // dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, - // spatial_dimension*spatial_dimension, "stress"); - // dumper.SetEmbeddedValue("displacements", 1); - // dumper.SetEmbeddedValue("applied_force", 1); - // dumper.SetEmbeddedValue("forces", 1); - // dumper.SetPrefix("paraview/"); - // dumper.Init(); - // dumper.Dump(); + iohelper::ElemType paraview_type = iohelper::TRIANGLE2; + UInt nb_element = mesh.getNbElement(type); + + /// initialize the paraview output + iohelper::DumperParaview dumper; + dumper.SetMode(iohelper::TEXT); + dumper.SetPoints(mesh.getNodes().values, + spatial_dimension, mesh.getNbNodes(), "explicit"); + dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, + paraview_type, nb_element, iohelper::C_MODE); + dumper.AddNodeDataField(model.getDisplacement().values, + spatial_dimension, "displacements"); + dumper.AddNodeDataField(model.getVelocity().values, + spatial_dimension, "velocity"); + dumper.AddNodeDataField(model.getAcceleration().values, + spatial_dimension, "acceleration"); + dumper.AddNodeDataField(model.getForce().values, + spatial_dimension, "applied_force"); + dumper.AddNodeDataField(model.getResidual().values, + spatial_dimension, "forces"); + dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, + spatial_dimension*spatial_dimension, "strain"); + dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, + spatial_dimension*spatial_dimension, "stress"); + dumper.SetEmbeddedValue("displacements", 1); + dumper.SetEmbeddedValue("applied_force", 1); + dumper.SetEmbeddedValue("forces", 1); + dumper.SetPrefix("paraview/"); + dumper.Init(); + dumper.Dump(); /// initial conditions Real loading_rate = 0.5; Real disp_update = loading_rate * time_step; for (UInt n = 0; n < nb_nodes; ++n) { velocity(n, 1) = loading_rate * position(n, 1); } // std::ofstream edis("edis.txt"); // std::ofstream erev("erev.txt"); // Vector & residual = model.getResidual(); // const Vector & stress = model.getMaterial(0).getStress(type); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { /// update displacement on extreme nodes for (UInt n = 0; n < nb_nodes; ++n) { if (position(n, 1) > 0.99 || position(n, 1) < -0.99) displacement(n, 1) += disp_update * position(n, 1); } model.checkCohesiveStress(); model.explicitPred(); model.updateResidual(); model.updateAcceleration(); model.explicitCorr(); if(s % 1 == 0) { - // dumper.SetPoints(mesh.getNodes().values, - // spatial_dimension, mesh.getNbNodes(), "explicit"); - // dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, - // paraview_type, nb_element, iohelper::C_MODE); - // dumper.AddNodeDataField(model.getDisplacement().values, - // spatial_dimension, "displacements"); - // dumper.AddNodeDataField(model.getVelocity().values, - // spatial_dimension, "velocity"); - // dumper.AddNodeDataField(model.getAcceleration().values, - // spatial_dimension, "acceleration"); - // dumper.AddNodeDataField(model.getForce().values, - // spatial_dimension, "applied_force"); - // dumper.AddNodeDataField(model.getResidual().values, - // spatial_dimension, "forces"); - // dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, - // spatial_dimension*spatial_dimension, "strain"); - // dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, - // spatial_dimension*spatial_dimension, "stress"); - // dumper.SetEmbeddedValue("displacements", 1); - // dumper.SetEmbeddedValue("applied_force", 1); - // dumper.SetEmbeddedValue("forces", 1); - // dumper.SetPrefix("paraview/"); - // dumper.Init(); - // dumper.Dump(); + dumper.SetPoints(mesh.getNodes().values, + spatial_dimension, mesh.getNbNodes(), "explicit"); + dumper.SetConnectivity((int *)mesh.getConnectivity(type).values, + paraview_type, nb_element, iohelper::C_MODE); + dumper.AddNodeDataField(model.getDisplacement().values, + spatial_dimension, "displacements"); + dumper.AddNodeDataField(model.getVelocity().values, + spatial_dimension, "velocity"); + dumper.AddNodeDataField(model.getAcceleration().values, + spatial_dimension, "acceleration"); + dumper.AddNodeDataField(model.getForce().values, + spatial_dimension, "applied_force"); + dumper.AddNodeDataField(model.getResidual().values, + spatial_dimension, "forces"); + dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, + spatial_dimension*spatial_dimension, "strain"); + dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, + spatial_dimension*spatial_dimension, "stress"); + dumper.SetEmbeddedValue("displacements", 1); + dumper.SetEmbeddedValue("applied_force", 1); + dumper.SetEmbeddedValue("forces", 1); + dumper.SetPrefix("paraview/"); + dumper.Init(); + dumper.Dump(); std::cout << "passing step " << s << "/" << max_steps << std::endl; } // Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); // Real Er = dynamic_cast (model.getMaterial(1)).getReversibleEnergy(); // edis << s << " " // << Ed << std::endl; // erev << s << " " // << Er << std::endl; } // edis.close(); // erev.close(); Real Ed = model.getEnergy("dissipated"); Real Edt = 200*std::sqrt(2); std::cout << Ed << " " << Edt << std::endl; if (Ed < Edt * 0.999 || Ed > Edt * 1.001) { std::cout << "The dissipated energy is incorrect" << std::endl; return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_extrinsic was passed!" << std::endl; return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/mesh.geo b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/mesh.geo index 4526cbf34..b7d030c25 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/mesh.geo +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/mesh.geo @@ -1,22 +1,13 @@ h = 1.0; -Point(1) = {0, 0, 0, h}; -Point(2) = {0, 1, 0, h}; -Point(3) = {-1, 1, 0, h}; -Point(4) = {-1, 0, 0, h}; -Point(6) = {1, 0, 0, h}; -Point(7) = {1, 1, 0, h}; -Point(8) = {1, -1, 0, h}; -Point(9) = {0, -1, 0, h}; -Point(10) = {-1, -1, 0, h}; -Line(1) = {10, 9}; -Line(2) = {9, 8}; -Line(3) = {8, 6}; -Line(4) = {6, 7}; -Line(5) = {7, 2}; -Line(6) = {2, 3}; -Line(7) = {3, 4}; -Line(8) = {4, 10}; -Line Loop(9) = {3, 4, 5, 6, 7, 8, 1, 2}; -Plane Surface(10) = {9}; -Physical Surface(11) = {10}; +Point(1) = {1, 1, 0, h}; +Point(2) = {-1, 1, 0, h}; +Point(3) = {-1, -1, 0, h}; +Point(4) = {1, -1, 0, h}; +Line(1) = {1, 2}; +Line(2) = {2, 3}; +Line(3) = {3, 4}; +Line(4) = {4, 1}; +Line Loop(5) = {4, 1, 2, 3}; +Plane Surface(6) = {5}; +Physical Surface(7) = {6}; diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc index 732aeb1eb..b986acf4f 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic/test_cohesive_intrinsic.cc @@ -1,259 +1,262 @@ /** * @file test_cohesive.cc * @author Marco Vocialta * @date Fri Feb 24 14:32:31 2012 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model_cohesive.hh" #include "material.hh" #include "io_helper.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; static void updateDisplacement(SolidMechanicsModelCohesive &, Vector &, ElementType, Real); int main(int argc, char *argv[]) { initialize(argc, argv); - // debug::setDebugLevel(dblDump); + debug::setDebugLevel(dblDump); const UInt spatial_dimension = 2; const UInt max_steps = 350; const ElementType type = _triangle_6; Mesh mesh(spatial_dimension); MeshIOMSH mesh_io; mesh_io.read("mesh.msh", mesh); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull("material.dat"); Real time_step = model.getStableTimeStep()*0.8; model.setTimeStep(time_step); // std::cout << "Time step: " << time_step << std::endl; model.assembleMassLumped(); /* ------------------------------------------------------------------------ */ /* Facet part */ /* ------------------------------------------------------------------------ */ - // std::cout << mesh << std::endl; + std::cout << mesh << std::endl; const Mesh & mesh_facets = model.getMeshFacets(); + std::cout << mesh_facets << std::endl; + const ElementType type_facet = mesh.getFacetElementType(type); UInt nb_facet = mesh_facets.getNbElement(type_facet); // const Vector & position = mesh.getNodes(); // Vector & displacement = model.getDisplacement(); // const Vector & connectivity = mesh_facets.getConnectivity(type_facet); Vector facet_insertion; Real * bary_facet = new Real[spatial_dimension]; for (UInt f = 0; f < nb_facet; ++f) { mesh_facets.getBarycenter(f, type_facet, bary_facet); - if (bary_facet[0] > -0.30 && bary_facet[0] < -0.20) facet_insertion.push_back(f); + if (bary_facet[0] > -0.26 && bary_facet[0] < -0.24) facet_insertion.push_back(f); } delete[] bary_facet; model.insertCohesiveElements(facet_insertion); // mesh_io.write("mesh_cohesive.msh", mesh); // std::cout << mesh << std::endl; /* ------------------------------------------------------------------------ */ /* End of facet part */ /* ------------------------------------------------------------------------ */ Vector & boundary = model.getBoundary(); // const Vector & residual = model.getResidual(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_element = mesh.getNbElement(type); /// boundary conditions for (UInt dim = 0; dim < spatial_dimension; ++dim) { for (UInt n = 0; n < nb_nodes; ++n) { boundary(n, dim) = true; } } model.updateResidual(); // iohelper::ElemType paraview_type = iohelper::TRIANGLE2; // /// initialize the paraview output // iohelper::DumperParaview dumper; // dumper.SetPoints(model.getFEM().getMesh().getNodes().values, // spatial_dimension, nb_nodes, "explicit"); // dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values, // paraview_type, nb_element, iohelper::C_MODE); // dumper.AddNodeDataField(model.getDisplacement().values, // spatial_dimension, "displacements"); // dumper.AddNodeDataField(model.getVelocity().values, // spatial_dimension, "velocity"); // dumper.AddNodeDataField(model.getAcceleration().values, // spatial_dimension, "acceleration"); // dumper.AddNodeDataField(model.getForce().values, // spatial_dimension, "applied_force"); // dumper.AddNodeDataField(model.getResidual().values, // spatial_dimension, "forces"); // dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, // spatial_dimension*spatial_dimension, "strain"); // dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, // spatial_dimension*spatial_dimension, "stress"); // dumper.SetEmbeddedValue("displacements", 1); // dumper.SetEmbeddedValue("applied_force", 1); // dumper.SetEmbeddedValue("forces", 1); // dumper.SetPrefix("paraview/"); // dumper.Init(); // dumper.Dump(); /// update displacement Vector elements; Real * bary = new Real[spatial_dimension]; for (UInt el = 0; el < nb_element; ++el) { mesh.getBarycenter(el, type, bary); if (bary[0] > -0.25) elements.push_back(el); } delete[] bary; Real increment = 0.01; updateDisplacement(model, elements, type, increment); // for (UInt n = 0; n < nb_nodes; ++n) { // if (position(n, 1) + displacement(n, 1) > 0) { // if (position(n, 0) == 0) { // displacement(n, 1) -= 0.25; // } // if (position(n, 0) == 1) { // displacement(n, 1) += 0.25; // } // } // } // std::ofstream edis("edis.txt"); // std::ofstream erev("erev.txt"); /// Main loop for (UInt s = 1; s <= max_steps; ++s) { model.explicitPred(); model.updateResidual(); model.updateAcceleration(); model.explicitCorr(); updateDisplacement(model, elements, type, increment); if(s % 1 == 0) { // dumper.Dump(); std::cout << "passing step " << s << "/" << max_steps << std::endl; } // // update displacement // for (UInt n = 0; n < nb_nodes; ++n) { // if (position(n, 1) + displacement(n, 1) > 0) { // displacement(n, 0) -= 0.01; // } // } // Real Ed = dynamic_cast (model.getMaterial(1)).getDissipatedEnergy(); // Real Er = dynamic_cast (model.getMaterial(1)).getReversibleEnergy(); // edis << s << " " // << Ed << std::endl; // erev << s << " " // << Er << std::endl; } // edis.close(); // erev.close(); Real Ed = model.getEnergy("dissipated"); Real Edt = 2 * sqrt(2); std::cout << Ed << " " << Edt << std::endl; if (Ed < Edt * 0.999 || Ed > Edt * 1.001) { std::cout << "The dissipated energy is incorrect" << std::endl; return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_intrinsic was passed!" << std::endl; return EXIT_SUCCESS; } static void updateDisplacement(SolidMechanicsModelCohesive & model, Vector & elements, ElementType type, Real increment) { Mesh & mesh = model.getFEM().getMesh(); UInt nb_element = elements.getSize(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); const Vector & connectivity = mesh.getConnectivity(type); Vector & displacement = model.getDisplacement(); Vector update(nb_nodes); + update.clear(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity(elements(el), n); if (!update(node)) { displacement(node, 0) += increment; // displacement(node, 1) += increment; update(node) = true; } } } } diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc index 68db41a71..c507eaca8 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc @@ -1,177 +1,180 @@ /** * @file test_interpolate_stress.cc * @author Marco Vocialta * @date Wed Jun 6 14:25:54 2012 * * @brief Test for the stress interpolation function * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "io_helper.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; static Real plane(Real, Real); int main(int argc, char *argv[]) { initialize(argc, argv); debug::setDebugLevel(dblWarning); const UInt spatial_dimension = 2; const ElementType type = _triangle_6; Mesh mesh(spatial_dimension); MeshIOMSH mesh_io; mesh_io.read("triangle.msh", mesh); const ElementType type_facet = mesh.getFacetElementType(type);; MeshUtils::buildAllFacets(mesh, mesh); SolidMechanicsModel model(mesh); /// model initialization model.initFull("../material.dat"); Vector & position = mesh.getNodes(); UInt nb_facet = mesh.getNbElement(type_facet); UInt nb_element = mesh.getNbElement(type); /// compute quadrature points positions on facets typedef FEMTemplate MyFEMType; model.registerFEMObject("FacetsFEM", mesh, spatial_dimension-1); model.getFEM("FacetsFEM").initShapeFunctions(); UInt nb_quad_per_facet = model.getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_tot_quad = nb_quad_per_facet * nb_facet; Vector quad_facets(nb_tot_quad, spatial_dimension); model.getFEM("FacetsFEM").interpolateOnQuadraturePoints(position, quad_facets, spatial_dimension, type_facet); Vector & facet_to_element = mesh.getSubelementToElement(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); - Vector el_q_facet(nb_element * nb_facet_per_elem * nb_quad_per_facet, - spatial_dimension); + ByElementTypeReal element_quad_facet; + element_quad_facet.alloc(nb_element * nb_facet_per_elem * nb_quad_per_facet, + spatial_dimension, + type); + + Vector & el_q_facet = element_quad_facet(type); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { UInt global_facet = facet_to_element(el, f).element; for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < spatial_dimension; ++s) { el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s) = quad_facets(global_facet * nb_quad_per_facet + q, s); } } } } /// compute quadrature points position of the elements UInt nb_quad_per_element = model.getFEM().getNbQuadraturePoints(type); UInt nb_tot_quad_el = nb_quad_per_element * nb_element; Vector quad_elements(nb_tot_quad_el, spatial_dimension); model.getFEM().interpolateOnQuadraturePoints(position, quad_elements, spatial_dimension, type); /// assign some values to stresses Vector & stress = const_cast&>(model.getMaterial(0).getStress(type)); for (UInt q = 0; q < nb_tot_quad_el; ++q) { for (UInt s = 0; s < spatial_dimension * spatial_dimension; ++s) { stress(q, s) = s * plane(quad_elements(q, 0), quad_elements(q, 1)); } } /// interpolate stresses on facets' quadrature points - model.getMaterial(0).initInterpolateElementalField(); + model.getMaterial(0).initElementalFieldInterpolation(element_quad_facet); Vector interpolated_stress(nb_element * nb_facet_per_elem * nb_quad_per_facet, stress.getNbComponent()); model.getMaterial(0).interpolateStress(type, - el_q_facet, interpolated_stress); /// check results for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < spatial_dimension * spatial_dimension; ++s) { Real x = el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, 0); Real y = el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, 1); Real theoretical = s * plane(x, y); Real numerical = interpolated_stress(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s); Real tolerance = 1000 * (theoretical + 1) * std::numeric_limits::epsilon(); if (std::abs(theoretical - numerical) > tolerance) { std::cout << "Theoretical and numerical values aren't coincident!" << std::endl; return EXIT_FAILURE; } } } } } std::cout << "OK: Stress interpolation test passed." << std::endl; return EXIT_SUCCESS; } Real plane(Real x, Real y) { return 1. + 2. * x + 3. * y; }