diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh index 96f4b5a0b..541c4855e 100644 --- a/src/common/aka_types.hh +++ b/src/common/aka_types.hh @@ -1,987 +1,995 @@ /** * @file aka_types.hh * * @author Nicolas Richart * * @date creation: Thu Feb 17 2011 * @date last modification: Tue Aug 19 2014 * * @brief description of the "simple" types * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_error.hh" #include "aka_fwd.hh" #include "aka_math.hh" #include "aka_array.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 enum NormType { L_1 = 1, L_2 = 2, L_inf = UInt(-1) }; template struct DimHelper { static inline void setDims(UInt m, UInt n, UInt p, UInt dims[dim]); }; template<> struct DimHelper<1> { static inline void setDims(UInt m, __attribute__((unused)) UInt n, __attribute__((unused)) UInt p, UInt dims[1]) { dims[0] = m; } }; template<> struct DimHelper<2> { static inline void setDims(UInt m, UInt n, __attribute__((unused)) UInt p, UInt dims[2]) { dims[0] = m; dims[1] = n; } }; template<> struct DimHelper<3> { static inline void setDims(UInt m, UInt n, UInt p, UInt dims[3]) { dims[0] = m; dims[1] = n; dims[2] = p; } }; /* -------------------------------------------------------------------------- */ template class TensorProxy { protected: TensorProxy(T * data, UInt m, UInt n, UInt p) { DimHelper::setDims(m, n, p, this->n); this->values = data; } TensorProxy(const TensorProxy & other) { this->values = other.storage(); for (UInt i = 0; i < ndim; ++i) this->n[i] = other.n[i]; } public: UInt size(UInt i) const { AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim << " dimensions, not " << (i + 1)); return n[i]; } inline UInt size() const{ UInt _size = 1; for (UInt d = 0; d < ndim; ++d) _size *= this->n[d]; return _size; } T * storage() const { return values; } // TensorProxy operator=(const TensorProxy & other) { // UInt _size = this->size(); // #ifndef AKANTU_NDEBUG // UInt _size_other = other.size(); // AKANTU_DEBUG_ASSERT(_size == _size_other, "The two tensor are not compatible in size"); // #endif // memcpy(this->values, other.storage(), _size * sizeof(T)); // return *this; // } protected: T * values; UInt n[ndim]; }; /* -------------------------------------------------------------------------- */ template class VectorProxy : public TensorProxy { typedef TensorProxy parent; public: VectorProxy(T * data = NULL, UInt n = 0) : parent(data, n, 0, 0) { } VectorProxy(const VectorProxy & src) : parent(src) { } // VectorProxy & operator=(const VectorProxy & other) { // parent::operator=(other); // return *this; // } }; template class MatrixProxy : public TensorProxy { typedef TensorProxy parent; public: MatrixProxy(T * data = NULL, UInt m = 0, UInt n = 0) : parent(data, m, n, 0) { } MatrixProxy(const MatrixProxy & src) : parent(src) { } // MatrixProxy & operator=(const MatrixProxy & other) { // parent::operator=(other); // return *this; // } }; template class Tensor3Proxy : public TensorProxy { typedef TensorProxy parent; public: Tensor3Proxy(T * data = NULL, UInt m = 0, UInt n = 0, UInt k = 0) : parent(data, m, n, k) { } Tensor3Proxy(const Tensor3Proxy & src) : parent(src) { } // Tensor3Proxy & operator=(const Tensor3Proxy & other) { // parent::operator=(other); // return *this; // } }; /* -------------------------------------------------------------------------- */ /* Tensor base class */ /* -------------------------------------------------------------------------- */ template class TensorStorage { public: typedef T value_type; protected: template void copySize(const TensorType & src) { for (UInt d = 0; d < ndim; ++d) this->n[d] = src.size(d); this->_size = src.size(); } TensorStorage() : values(NULL), wrapped(false) { for (UInt d = 0; d < ndim; ++d) this->n[d] = 0; _size = 0; } TensorStorage(const TensorProxy & proxy) { this->copySize(proxy); this->values = proxy.storage(); this->wrapped = true; } protected: TensorStorage(const TensorStorage & src) { } public: TensorStorage(const TensorStorage & src, bool deep_copy) : values(NULL), wrapped(false) { if(deep_copy) this->deepCopy(src); else this->shallowCopy(src); } protected: TensorStorage(UInt m, UInt n, UInt p, const T & def) { DimHelper::setDims(m, n, p, this->n); this->computeSize(); this->values = new T[this->_size]; this->set(def); this->wrapped = false; } TensorStorage(T * data, UInt m, UInt n, UInt p) { DimHelper::setDims(m, n, p, this->n); this->computeSize(); this->values = data; this->wrapped = true; } public: /* ------------------------------------------------------------------------ */ template inline void shallowCopy(const TensorType & src) { this->copySize(src); if(!this->wrapped) delete[] this->values; this->values = src.storage(); this->wrapped = true; } /* ------------------------------------------------------------------------ */ template inline void deepCopy(const TensorType & src) { this->copySize(src); if(!this->wrapped) delete [] this->values; this->values = new T[this->_size]; memcpy(this->values, src.storage(), this->_size * sizeof(T)); this->wrapped = false; } virtual ~TensorStorage() { if(!this->wrapped) delete [] this->values; } inline TensorStorage & operator=(const RetType & src) { if(this != &src) { if (this->wrapped) { AKANTU_DEBUG_ASSERT(this->_size == src.size(), "vectors of different size"); memcpy(this->values, src.storage(), this->_size * sizeof(T)); } else { deepCopy(src); } } return *this; } /* ------------------------------------------------------------------------ */ template inline RetType & operator+=(const TensorStorage & other) { T * a = this->storage(); T * b = other.storage(); AKANTU_DEBUG_ASSERT(_size == other.size(), "The two tensors do not have the same size, they cannot be subtracted"); for (UInt i = 0; i < _size; ++i) *(a++) += *(b++); return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ template inline RetType & operator-=(const TensorStorage & other) { T * a = this->storage(); T * b = other.storage(); AKANTU_DEBUG_ASSERT(_size == other.size(), "The two tensors do not have the same size, they cannot be subtracted"); for (UInt i = 0; i < _size; ++i) *(a++) -= *(b++); return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ inline RetType & operator+=(const T & x) { T * a = this->values; for (UInt i = 0; i < _size; ++i) *(a++) += x; return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ inline RetType & operator-=(const T & x) { T * a = this->values; for (UInt i = 0; i < _size; ++i) *(a++) -= x; return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ inline RetType & operator*=(const T & x) { T * a = this->storage(); for (UInt i = 0; i < _size; ++i) *(a++) *= x; return *(static_cast(this)); } /* ---------------------------------------------------------------------- */ inline RetType & operator/=(const T & x) { T * a = this->values; for (UInt i = 0; i < _size; ++i) *(a++) /= x; return *(static_cast(this)); } /* ------------------------------------------------------------------------ */ T * storage() const { return values; } UInt size() const { return _size; } UInt size(UInt i) const { AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim << " dimensions, not " << (i + 1)); return n[i]; }; /* ------------------------------------------------------------------------ */ inline void clear() { memset(values, 0, _size * sizeof(T)); }; inline void set(const T & t) { std::fill_n(values, _size, t); }; template inline void copy(const TensorType & other) { AKANTU_DEBUG_ASSERT(_size == other.size(), "The two tensors do not have the same size, they cannot be copied"); memcpy(values, other.storage(), _size * sizeof(T)); } protected: friend class Array; inline void computeSize() { _size = 1; for (UInt d = 0; d < ndim; ++d) _size *= this->n[d]; } protected: template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm += std::pow(std::abs(*it), norm_type); return std::pow(_norm, 1./norm_type); } }; template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm += std::abs(*it); return _norm; } }; template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm += *it * *it; return sqrt(_norm); } }; template struct NormHelper { template static R norm(const Ten & ten) { R _norm = 0.; R * it = ten.storage(); R * end = ten.storage() + ten.size(); for (; it < end; ++it) _norm = std::max(std::abs(*it), _norm); return _norm; } }; public: /*----------------------------------------------------------------------- */ /// "Entrywise" norm norm @f[ \|\boldsymbol{T}\|_p = \left( /// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}} /// @f] template inline T norm() const { return NormHelper::norm(*this); } protected: UInt n[ndim]; UInt _size; T * values; bool wrapped; }; /* -------------------------------------------------------------------------- */ /* Vector */ /* -------------------------------------------------------------------------- */ template class Vector : public TensorStorage< T, 1, Vector > { typedef TensorStorage< T, 1, Vector > parent; public: typedef typename parent::value_type value_type; public: Vector() : parent() {} Vector(UInt n, const T & def = T()) : parent(n, 0, 0, def) { } Vector(T * data, UInt n) : parent(data, n, 0, 0) { } Vector(const Vector & src, bool deep_copy = true) : parent(src, deep_copy) { } Vector(const VectorProxy & src) : parent(src) { } public: virtual ~Vector() { }; /* ------------------------------------------------------------------------ */ inline Vector & operator=(const Vector & src) { parent::operator=(src); return *this; } /* ------------------------------------------------------------------------ */ inline T& operator()(UInt i) { return *(this->values + i); }; inline const T& operator()(UInt i) const { return *(this->values + i); }; inline T& operator[](UInt i) { return *(this->values + i); }; inline const T& operator[](UInt i) const { return *(this->values + i); }; /* ------------------------------------------------------------------------ */ inline Vector & operator*=(Real x) { return parent::operator*=(x); } inline Vector & operator/=(Real x) { return parent::operator/=(x); } /* ------------------------------------------------------------------------ */ inline Vector & operator*=(const Vector & vect) { T * a = this->storage(); T * b = vect.storage(); for (UInt i = 0; i < this->_size; ++i) *(a++) *= *(b++); return *this; } /* ------------------------------------------------------------------------ */ inline Real dot(const Vector & vect) const { return Math::vectorDot(this->values, vect.storage(), this->_size); } + /* ------------------------------------------------------------------------ */ + inline Real mean() const { + Real mean = 0; + T * a = this->storage(); + for (UInt i = 0; i < this->_size; ++i) mean = *(a++); + return mean / this->_size; + } + /* ------------------------------------------------------------------------ */ inline Vector & crossProduct(const Vector & v1, const Vector & v2) { AKANTU_DEBUG_ASSERT(this->size() == 3, "crossProduct is only defined in 3D (n=" << this->size() << ")"); AKANTU_DEBUG_ASSERT(this->size() == v1.size() && this->size() == v2.size(), "crossProduct is not a valid operation non matching size vectors"); Math::vectorProduct3(v1.storage(), v2.storage(), this->values); return *this; } /* ------------------------------------------------------------------------ */ inline void solve(Matrix & A, const Vector & b) { AKANTU_DEBUG_ASSERT(this->size() == A.rows() && this->_size = A.cols(), "The solution vector as a mismatch in size with the matrix"); AKANTU_DEBUG_ASSERT(this->_size == b._size, "The rhs vector as a mismatch in size with the matrix"); Math::solve(this->_size, A.storage(), this->values, b.storage()); } /* ------------------------------------------------------------------------ */ template inline void mul(const Matrix & A, const Vector & x, Real alpha = 1.0); /* ------------------------------------------------------------------------ */ inline Real norm() const { return parent::template norm(); } template inline Real norm() const { return parent::template norm(); } /* ------------------------------------------------------------------------ */ inline void normalize() { Real n = norm(); operator/=(n); } /* ------------------------------------------------------------------------ */ /// norm of (*this - x) inline Real distance(const Vector & y) const { Real * vx = this->values; Real * vy = y.storage(); Real sum_2 = 0; for (UInt i = 0; i < this->_size; ++i, ++vx, ++vy) sum_2 += (*vx - *vy)*(*vx - *vy); return sqrt(sum_2); } /* ------------------------------------------------------------------------ */ inline bool equal(const Vector & v, Real tolerance = Math::getTolerance()) const { T * a = this->storage(); T * b = v.storage(); UInt i = 0; while (i < this->_size && (std::abs(*(a++) - *(b++)) < tolerance)) ++i; return i == this->_size; } /* ------------------------------------------------------------------------ */ inline short compare(const Vector & v, Real tolerance = Math::getTolerance()) const { T * a = this->storage(); T * b = v.storage(); for (UInt i(0); i < this->_size; ++i, ++a, ++b) { if(std::abs(*a - *b) > tolerance) return (((*a - *b) > tolerance) ? 1 : -1); } return 0; } /* ------------------------------------------------------------------------ */ inline bool operator==(const Vector & v) const { return equal(v); } inline bool operator<(const Vector & v) const { return compare(v) == -1; } inline bool operator>(const Vector & v) const { return compare(v) == 1; } /* ------------------------------------------------------------------------ */ /// 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 << "Vector<" << debug::demangle(typeid(T).name()) << ">(" << this->_size <<") : ["; for (UInt i = 0; i < this->_size; ++i) { if(i != 0) stream << ", "; stream << this->values[i]; } stream << "]"; } friend class ::akantu::Array; }; 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.size()); r = a; r *= scalar; return r; } template Vector operator+(const Vector & a, const Vector & b) { Vector r(a.size()); r = a; r += b; return r; } template Vector operator-(const Vector& a, const Vector& b) { Vector r(a.size()); r = a; r -= b; return r; } template Matrix operator*(T scalar, const Matrix& a) { Matrix r(a.rows(), a.cols()); r = a; r *= scalar; return r; } template Matrix operator+(const Matrix& a, const Matrix& b) { Matrix r(a.rows(), a.cols()); r = a; r += b; return r; } template Matrix operator-(const Matrix& a, const Matrix& b) { Matrix r(a.rows(), a.cols()); r = a; r -= b; return r; } /* ------------------------------------------------------------------------ */ template<> inline bool Vector::equal(const Vector & v, __attribute__((unused)) Real tolerance) const { UInt * a = this->storage(); UInt * b = v.storage(); UInt i = 0; while (i < this->_size && (*(a++) == *(b++))) ++i; return i == this->_size; } /* ------------------------------------------------------------------------ */ /* Matrix */ /* ------------------------------------------------------------------------ */ template class Matrix : public TensorStorage< T, 2, Matrix > { typedef TensorStorage< T, 2, Matrix > parent; public: typedef typename parent::value_type value_type; public: Matrix() : parent() {} Matrix(UInt m, UInt n, const T & def = T()) : parent(m, n, 0, def) { } Matrix(T * data, UInt m, UInt n) : parent(data, m, n, 0) { } Matrix(const Matrix & src, bool deep_copy = true) : parent(src, deep_copy) { } Matrix(const MatrixProxy & src) : parent(src) { } virtual ~Matrix() { } /* ------------------------------------------------------------------------ */ inline Matrix & operator=(const Matrix & src) { parent::operator=(src); return *this; } public: /* ---------------------------------------------------------------------- */ UInt rows() const { return this->n[0]; } UInt cols() const { return this->n[1]; } /* ---------------------------------------------------------------------- */ inline T& operator()(UInt i, UInt j) { return this->at(i,j); } inline T& at(UInt i, UInt j) { return *(this->values + i + j*this->n[0]); } inline const T& operator()(UInt i, UInt j) const { return this->at(i,j); } inline T& at(UInt i, UInt j) const { return *(this->values + i + j*this->n[0]); } /// give a line vector wrapped on the column i inline VectorProxy operator()(UInt j) { return VectorProxy(this->values + j*this->n[0], this->n[0]); } inline const VectorProxy operator()(UInt j) const { return VectorProxy(this->values + j*this->n[0], this->n[0]); } inline T& operator[](UInt idx) { return *(this->values + idx); }; inline const T& operator[](UInt idx) const { return *(this->values + idx); }; /* ---------------------------------------------------------------------- */ inline Matrix operator* (const Matrix & B) { Matrix C(this->rows(), B.cols()); C.mul(*this, B); return C; } /* ----------------------------------------------------------------------- */ inline Matrix & operator*= (const T & x) { return parent::operator*= (x); } inline Matrix & operator*= (const Matrix & B) { Matrix C(*this); this->mul(C, B); return *this; } /* ---------------------------------------------------------------------- */ template inline void mul(const Matrix & A, const Matrix & B, T alpha = 1.0) { UInt k = A.cols(); if(tr_A) k = A.rows(); #ifndef AKANTU_NDEBUG if (tr_B){ AKANTU_DEBUG_ASSERT(k == B.cols(), "matrices to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->cols() == B.rows(), "matrices to multiply have no fit dimensions"); } else { AKANTU_DEBUG_ASSERT(k == B.rows(), "matrices to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->cols() == B.cols(), "matrices to multiply have no fit dimensions"); } if (tr_A){ AKANTU_DEBUG_ASSERT(this->rows() == A.cols(), "matrices to multiply have no fit dimensions"); } else{ AKANTU_DEBUG_ASSERT(this->rows() == A.rows(), "matrices to multiply have no fit dimensions"); } #endif //AKANTU_NDEBUG Math::matMul(this->rows(), this->cols(), k, alpha, A.storage(), B.storage(), 0., this->storage()); } /* ---------------------------------------------------------------------- */ inline void outerProduct(const Vector & A, const Vector & B) { AKANTU_DEBUG_ASSERT(A.size() == this->rows() && B.size() == this->cols(), "A and B are not compatible with the size of the matrix"); for (UInt i = 0; i < this->rows(); ++i) { for (UInt j = 0; j < this->cols(); ++j) { this->values[i + j * this->rows()] += A[i] * B[j]; } } } private: class EigenSorter { public: EigenSorter(const Vector & eigs) : eigs(eigs) { } bool operator()(const UInt & a, const UInt & b) const { return (eigs(a) > eigs(b)); } private: const Vector & eigs; }; public: /* ---------------------------------------------------------------------- */ inline void eig(Vector & eigenvalues, Matrix & eigenvectors) const { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "eig is not a valid operation on a rectangular matrix"); AKANTU_DEBUG_ASSERT(eigenvalues.size() == this->cols(), "eigenvalues should be of size " << this->cols() << "."); #ifndef AKANTU_NDEBUG if(eigenvectors.storage() != NULL) AKANTU_DEBUG_ASSERT((eigenvectors.rows() == eigenvectors.cols()) && (eigenvectors.rows() == this->cols()), "Eigenvectors needs to be a square matrix of size " << this->cols() << " x " << this->cols() << "."); #endif Matrix tmp = *this; Vector tmp_eigs(eigenvalues.size()); Matrix tmp_eig_vects(eigenvectors.rows(), eigenvectors.cols()); if(tmp_eig_vects.rows() == 0 || tmp_eig_vects.cols() == 0) Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage()); else Math::matrixEig(tmp.cols(), tmp.storage(), tmp_eigs.storage(), tmp_eig_vects.storage()); Vector perm(eigenvalues.size()); for (UInt i = 0; i < perm.size(); ++i) perm(i) = i; std::sort(perm.storage(), perm.storage() + perm.size(), EigenSorter(tmp_eigs)); for (UInt i = 0; i < perm.size(); ++i) eigenvalues(i) = tmp_eigs(perm(i)); if(tmp_eig_vects.rows() != 0 && tmp_eig_vects.cols() != 0) for (UInt i = 0; i < perm.size(); ++i) { for (UInt j = 0; j < eigenvectors.rows(); ++j) { eigenvectors(j, i) = tmp_eig_vects(j, perm(i)); } } } /* ---------------------------------------------------------------------- */ inline void eig(Vector & eigenvalues) const { Matrix empty; eig(eigenvalues, empty); } /* ---------------------------------------------------------------------- */ inline void eye(T alpha = 1.) { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "eye is not a valid operation on a rectangular matrix"); this->clear(); for (UInt i = 0; i < this->cols(); ++i) { this->values[i + i * this->rows()] = alpha; } } /* ---------------------------------------------------------------------- */ static inline Matrix eye(UInt m, T alpha = 1.) { Matrix tmp(m, m); tmp.eye(alpha); return tmp; } /* ---------------------------------------------------------------------- */ inline T trace() const { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "trace is not a valid operation on a rectangular matrix"); T trace = 0.; for (UInt i = 0; i < this->rows(); ++i) { trace += this->values[i + i * this->rows()]; } return trace; } /* ---------------------------------------------------------------------- */ inline Matrix transpose() const { Matrix tmp(this->cols(), this->rows()); for (UInt i = 0; i < this->rows(); ++i) { for (UInt j = 0; j < this->cols(); ++j) { tmp(j,i) = operator()(i, j); } } return tmp; } /* ---------------------------------------------------------------------- */ inline void inverse(const Matrix & A) { AKANTU_DEBUG_ASSERT(A.cols() == A.rows(), "inv is not a valid operation on a rectangular matrix"); AKANTU_DEBUG_ASSERT(this->cols() == A.cols(), "the matrix should have the same size as its inverse"); if(this->cols() == 1) *this->values = 1./ *A.storage(); else if(this->cols() == 2) Math::inv2(A.storage(), this->values); else if(this->cols() == 3) Math::inv3(A.storage(), this->values); else Math::inv(this->cols(), A.storage(), this->values); } /* --------------------------------------------------------------------- */ inline T det() const { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "inv is not a valid operation on a rectangular matrix"); if(this->cols() == 1) return *(this->values); else if(this->cols() == 2) return Math::det2(this->values); else if(this->cols() == 3) return Math::det3(this->values); else return Math::det(this->cols(), this->values); } /* --------------------------------------------------------------------- */ inline T doubleDot(const Matrix & other) const { AKANTU_DEBUG_ASSERT(this->cols() == this->rows(), "doubleDot is not a valid operation on a rectangular matrix"); if(this->cols() == 1) return *(this->values) * *(other.storage()); else if(this->cols() == 2) return Math::matrixDoubleDot22(this->values, other.storage()); else if(this->cols() == 3) return Math::matrixDoubleDot33(this->values, other.storage()); else AKANTU_DEBUG_ERROR("doubleDot is not defined for other spatial dimensions" << " than 1, 2 or 3."); return T(); } /* ---------------------------------------------------------------------- */ /// 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 << "Matrix<" << debug::demangle(typeid(T).name()) << ">(" << this->n[0] << "," << this->n[1] <<") :" << "["; for (UInt i = 0; i < this->n[0]; ++i) { if(i != 0) stream << ", "; stream << "["; for (UInt j = 0; j < this->n[1]; ++j) { if(j != 0) stream << ", "; stream << operator()(i, j); } stream << "]"; } stream << "]"; }; }; /* ------------------------------------------------------------------------ */ template template inline void Vector::mul(const Matrix & A, const Vector & x, Real alpha) { #ifndef AKANTU_NDEBUG UInt n = x.size(); if (tr_A){ AKANTU_DEBUG_ASSERT(n == A.rows(), "matrix and vector to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->size() == A.cols(), "matrix and vector to multiply have no fit dimensions"); } else { AKANTU_DEBUG_ASSERT(n == A.cols(), "matrix and vector to multiply have no fit dimensions"); AKANTU_DEBUG_ASSERT(this->size() == A.rows(), "matrix and vector to multiply have no fit dimensions"); } #endif Math::matVectMul(A.rows(), A.cols(), alpha, A.storage(), x.storage(), 0., this->storage()); } /* -------------------------------------------------------------------------- */ template 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; } /* ------------------------------------------------------------------------ */ /* Tensor3 */ /* ------------------------------------------------------------------------ */ template class Tensor3 : public TensorStorage< T, 3, Tensor3 > { typedef TensorStorage< T, 3, Tensor3 > parent; public: typedef typename parent::value_type value_type; public: Tensor3() : parent() { }; Tensor3(UInt m, UInt n, UInt p, const T & def = T()) : parent(m, n, p, def) { } Tensor3(T * data, UInt m, UInt n, UInt p) : parent(data, m, n, p) { } Tensor3(const Tensor3 & src, bool deep_copy = true) : parent(src, deep_copy) { } public: /* ------------------------------------------------------------------------ */ inline Tensor3 & operator=(const Tensor3 & src) { parent::operator=(src); return *this; } /* ---------------------------------------------------------------------- */ inline T& operator()(UInt i, UInt j, UInt k) { return *(this->values + (k*this->n[0] + i)*this->n[1] + j); }; inline const T& operator()(UInt i, UInt j, UInt k) const { return *(this->values + (k*this->n[0] + i)*this->n[1] + j); }; inline MatrixProxy operator()(UInt k) { return MatrixProxy(this->values + k*this->n[0]*this->n[1], this->n[0], this->n[1]); } inline const MatrixProxy operator()(UInt k) const { return MatrixProxy(this->values + k*this->n[0]*this->n[1], this->n[0], this->n[1]); } inline MatrixProxy operator[](UInt k) { return Matrix(this->values + k*this->n[0]*this->n[1], this->n[0], this->n[1]); } inline const MatrixProxy operator[](UInt k) const { return MatrixProxy(this->values + k*this->n[0]*this->n[1], this->n[0], this->n[1]); } }; __END_AKANTU__ #endif /* __AKANTU_AKA_TYPES_HH__ */ diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index 2d4935129..8a450aaf8 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,2033 +1,2049 @@ /** * @file mesh_utils.cc * * @author Guillaume Anciaux * @author Leonardo Snozzi * @author Marco Vocialta * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * * @date creation: Fri Aug 20 2010 * @date last modification: Mon Jun 09 2014 * * @brief All mesh utils necessary for various tasks * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_utils.hh" #include "aka_safe_enum.hh" #include "fe_engine.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == _all_dimensions) spatial_dimension = mesh.getSpatialDimension(); /// count number of occurrence of each node UInt nb_nodes = mesh.getNbNodes(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); AKANTU_DEBUG_ASSERT(mesh.firstType(spatial_dimension) != mesh.lastType(spatial_dimension), "Some elements must be found in right dimension to compute facets!"); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt, _ek_not_defined); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_not_defined); for (; first != last; ++first) { ElementType type = *first; UInt nb_element = mesh.getNbElement(type, *gt); Array::const_iterator< Vector > conn_it = mesh.getConnectivity(type, *gt).begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) for (UInt n = 0; n < conn_it->size(); ++n) ++node_to_elem.rowOffset((*conn_it)(n)); } } node_to_elem.countToCSR(); node_to_elem.resizeCols(); /// rearrange element to get the node-element list Element e; node_to_elem.beginInsertions(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt, _ek_not_defined); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_not_defined); e.ghost_type = *gt; for (; first != last; ++first) { ElementType type = *first; e.type = type; e.kind = Mesh::getKind(type); UInt nb_element = mesh.getNbElement(type, *gt); Array::const_iterator< Vector > conn_it = mesh.getConnectivity(type, *gt).begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; for (UInt n = 0; n < conn_it->size(); ++n) node_to_elem.insertInRow((*conn_it)(n), e); } } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * This function should disappear in the future (used in mesh partitioning) */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == _all_dimensions) spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes = mesh.getNbNodes(); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); UInt nb_good_types = 0; UInt nb_nodes_per_element[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types]; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type); conn_val[nb_good_types] = mesh.getConnectivity(type, _not_ghost).storage(); nb_element[nb_good_types] = mesh.getConnectivity(type, _not_ghost).getSize(); nb_good_types++; } AKANTU_DEBUG_ASSERT(nb_good_types != 0, "Some elements must be found in right dimension to compute facets!"); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); /// count number of occurrence of each node for (UInt t = 0; t < nb_good_types; ++t) { for (UInt el = 0; el < nb_element[t]; ++el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) { ++node_to_elem.rowOffset(conn_val[t][el_offset + n]); } } } node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// rearrange element to get the node-element list for (UInt t = 0, linearized_el = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linearized_el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) node_to_elem.insertInRow(conn_val[t][el_offset + n], linearized_el); } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsElementTypeMap(const Mesh & mesh, CSR & node_to_elem, const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_elements = mesh.getConnectivity(type, ghost_type).getSize(); UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage(); /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); /// count number of occurrence of each node for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el*nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) ++node_to_elem.rowOffset(conn_val[el_offset + n]); } /// convert the occurrence array in a csr one node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// save the element index in the node-element list for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el*nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { node_to_elem.insertInRow(conn_val[el_offset + n], el); } } node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacets(Mesh & mesh){ AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { mesh.getConnectivity(*it, *gt).resize(0); // \todo inform the mesh event handler } } buildFacetsDimension(mesh, mesh, true, spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt to_dimension, DistributedSynchronizer * synchronizer) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); buildAllFacets(mesh, mesh_facets, spatial_dimension, to_dimension, synchronizer); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(const Mesh & mesh, Mesh & mesh_facets, UInt from_dimension, UInt to_dimension, DistributedSynchronizer * synchronizer) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mesh_facets.isMeshFacets(), "The mesh_facets should be initialized with initMeshFacets"); const ElementTypeMapArray * prank_to_element = NULL; if (synchronizer) { synchronizer->buildPrankToElement(); prank_to_element = &synchronizer->getPrankToElement(); } /// generate facets buildFacetsDimension(mesh, mesh_facets, false, from_dimension, prank_to_element); /// copy nodes type mesh_facets.nodes_type.resize(mesh.nodes_type.getSize()); mesh_facets.nodes_type.copy(mesh.nodes_type); /// sort facets and generate subfacets for (UInt i = from_dimension - 1; i > to_dimension; --i) { buildFacetsDimension(mesh_facets, mesh_facets, false, i, prank_to_element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacetsDimension(const Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension, const ElementTypeMapArray * prank_to_element){ AKANTU_DEBUG_IN(); // save the current parent of mesh_facets and set it temporarly to mesh since // mesh is the one containing the elements for which mesh_facets has the subelements // example: if the function is called with mesh = mesh_facets const Mesh & mesh_facets_parent = mesh_facets.getMeshParent(); mesh_facets.defineMeshParent(mesh); UInt spatial_dimension = mesh.getSpatialDimension(); const Array & mesh_facets_nodes = mesh_facets.getNodes(); const Array::const_vector_iterator mesh_facets_nodes_it = mesh_facets_nodes.begin(spatial_dimension); CSR node_to_elem; buildNode2Elements(mesh, node_to_elem, dimension); Array counter; + std::vector connected_elements; + + // init the SubelementToElement data to improve performance + for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { + GhostType ghost_type = *gt; + Mesh::type_iterator first = mesh.firstType(dimension, ghost_type); + Mesh::type_iterator last = mesh.lastType(dimension, ghost_type); + + for(; first != last; ++first) { + ElementType type = *first; + mesh_facets.getSubelementToElementPointer(type, ghost_type); + } + } Element current_element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; GhostType facet_ghost_type = ghost_type; current_element.ghost_type = ghost_type; Mesh::type_iterator first = mesh.firstType(dimension, ghost_type); Mesh::type_iterator last = mesh.lastType(dimension, ghost_type); for(; first != last; ++first) { ElementType type = *first; ElementType facet_type = mesh.getFacetType(type); current_element.type = type; UInt nb_element = mesh.getNbElement(type, ghost_type); Array< std::vector > * element_to_subelement = mesh_facets.getElementToSubelementPointer(facet_type, ghost_type); Array * connectivity_facets = mesh_facets.getConnectivityPointer(facet_type, ghost_type); + UInt nb_facet_per_element = mesh.getNbFacetsPerElement(type); + const Array & element_connectivity = mesh.getConnectivity(type, ghost_type); + const Matrix facet_local_connectivity = mesh.getFacetLocalConnectivity(type); + UInt nb_nodes_per_facet = connectivity_facets->getNbComponent(); + Vector facet(nb_nodes_per_facet); for (UInt el = 0; el < nb_element; ++el) { current_element.element = el; - Matrix facets = mesh.getFacetConnectivity(el, type, ghost_type); - UInt nb_nodes_per_facet = facets.cols(); - for (UInt f = 0; f < facets.rows(); ++f) { - Vector facet(nb_nodes_per_facet); - for (UInt n = 0; n < nb_nodes_per_facet; ++n) facet(n) = facets(f, n); + for (UInt f = 0; f < nb_facet_per_element; ++f) { + for (UInt n = 0; n < nb_nodes_per_facet; ++n) + facet(n) = element_connectivity(el, facet_local_connectivity(f, n)); - UInt first_node_nb_elements = node_to_elem.getNbCols(facets(f, 0)); + UInt first_node_nb_elements = node_to_elem.getNbCols(facet(0)); counter.resize(first_node_nb_elements); counter.clear(); //loop over the other nodes to search intersecting elements, //which are the elements that share another node with the //starting element after first_node CSR::iterator first_node_elements = node_to_elem.begin(facet(0)); CSR::iterator first_node_elements_end = node_to_elem.end(facet(0)); UInt local_el = 0; for (; first_node_elements != first_node_elements_end; ++first_node_elements, ++local_el) { for (UInt n = 1; n < nb_nodes_per_facet; ++n) { CSR::iterator node_elements_begin = node_to_elem.begin(facet(n)); CSR::iterator node_elements_end = node_to_elem.end (facet(n)); counter(local_el) += std::count(node_elements_begin, node_elements_end, *first_node_elements); } } // counting the number of elements connected to the facets and // taking the minimum element number, because the facet should // be inserted just once UInt nb_element_connected_to_facet = 0; Element minimum_el = ElementNull; - Array connected_elements; + connected_elements.resize(0); for (UInt el_f = 0; el_f < first_node_nb_elements; el_f++) { Element real_el = node_to_elem(facet(0), el_f); if (counter(el_f) == nb_nodes_per_facet - 1) { ++nb_element_connected_to_facet; minimum_el = std::min(minimum_el, real_el); connected_elements.push_back(real_el); } } if (minimum_el == current_element) { bool full_ghost_facet = false; UInt n = 0; while (n < nb_nodes_per_facet && mesh.isPureGhostNode(facet(n))) ++n; if (n == nb_nodes_per_facet) full_ghost_facet = true; if (!full_ghost_facet) { if (!boundary_only || (boundary_only && nb_element_connected_to_facet == 1)) { std::vector elements; // build elements_on_facets: linearized_el must come first // in order to store the facet in the correct direction // and avoid to invert the sign in the normal computation elements.push_back(current_element); /// boundary facet if (nb_element_connected_to_facet == 1) elements.push_back(ElementNull); /// internal facet else if (nb_element_connected_to_facet == 2) { - elements.push_back(connected_elements(1)); + elements.push_back(connected_elements[1]); /// check if facet is in between ghost and normal /// elements: if it's the case, the facet is either /// ghost or not ghost. The criterion to decide this /// is arbitrary. It was chosen to check the processor /// id (prank) of the two neighboring elements. If /// prank of the ghost element is lower than prank of /// the normal one, the facet is not ghost, otherwise /// it's ghost GhostType gt[2] = { _not_ghost, _not_ghost }; - for (UInt el = 0; el < connected_elements.getSize(); ++el) - gt[el] = connected_elements(el).ghost_type; + for (UInt el = 0; el < connected_elements.size(); ++el) + gt[el] = connected_elements[el].ghost_type; if (gt[0] + gt[1] == 1) { if (prank_to_element) { UInt prank[2]; for (UInt el = 0; el < 2; ++el) { - UInt current_el = connected_elements(el).element; - ElementType current_type = connected_elements(el).type; - GhostType current_gt = connected_elements(el).ghost_type; + UInt current_el = connected_elements[el].element; + ElementType current_type = connected_elements[el].type; + GhostType current_gt = connected_elements[el].ghost_type; const Array & prank_to_el = (*prank_to_element)(current_type, current_gt); prank[el] = prank_to_el(current_el); } bool ghost_one = (gt[0] != _ghost); if (prank[ghost_one] > prank[!ghost_one]) facet_ghost_type = _not_ghost; else facet_ghost_type = _ghost; connectivity_facets = mesh_facets.getConnectivityPointer(facet_type, facet_ghost_type); element_to_subelement = mesh_facets.getElementToSubelementPointer(facet_type, facet_ghost_type); } } } /// facet of facet else { for (UInt i = 1; i < nb_element_connected_to_facet; ++i) { - elements.push_back(connected_elements(i)); + elements.push_back(connected_elements[i]); } } element_to_subelement->push_back(elements); connectivity_facets->push_back(facet); /// current facet index UInt current_facet = connectivity_facets->getSize() - 1; /// loop on every element connected to current facet and /// insert current facet in the first free spot of the /// subelement_to_element vector for (UInt elem = 0; elem < elements.size(); ++elem) { Element loc_el = elements[elem]; if (loc_el.type != _not_defined) { Array & subelement_to_element = - *mesh_facets.getSubelementToElementPointer(type, loc_el.ghost_type); + mesh_facets.getSubelementToElement(type, loc_el.ghost_type); - for (UInt f_in = 0; f_in < facets.rows(); ++f_in) { + for (UInt f_in = 0; f_in < nb_facet_per_element; ++f_in) { if (subelement_to_element(loc_el.element, f_in).type == _not_defined) { subelement_to_element(loc_el.element, f_in).type = facet_type; subelement_to_element(loc_el.element, f_in).element = current_facet; subelement_to_element(loc_el.element, f_in).ghost_type = facet_ghost_type; break; } } } } /// reset connectivity in case a facet was found in /// between ghost and normal elements if (facet_ghost_type != ghost_type) { facet_ghost_type = ghost_type; connectivity_facets = mesh_facets.getConnectivityPointer(facet_type, facet_ghost_type); element_to_subelement = mesh_facets.getElementToSubelementPointer(facet_type, facet_ghost_type); } } } } } } } } // restore the parent of mesh_facet mesh_facets.defineMeshParent(mesh_facets_parent); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, UInt * local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Array & old_nodes_numbers) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::map renumbering_map; for (UInt i = 0; i < old_nodes_numbers.getSize(); ++i) { renumbering_map[old_nodes_numbers(i)] = i; } /// renumber the nodes renumberNodesInConnectivity(local_connectivities, (nb_local_element + nb_ghost_element)*nb_nodes_per_element, renumbering_map); std::map::iterator it = renumbering_map.begin(); std::map::iterator end = renumbering_map.end(); old_nodes_numbers.resize(renumbering_map.size()); for (;it != end; ++it) { old_nodes_numbers(it->second) = it->first; } renumbering_map.clear(); /// copy the renumbered connectivity to the right place Array * local_conn = mesh.getConnectivityPointer(type); local_conn->resize(nb_local_element); memcpy(local_conn->storage(), local_connectivities, nb_local_element * nb_nodes_per_element * sizeof(UInt)); Array * ghost_conn = mesh.getConnectivityPointer(type, _ghost); ghost_conn->resize(nb_ghost_element); memcpy(ghost_conn->storage(), local_connectivities + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberNodesInConnectivity(UInt * list_nodes, UInt nb_nodes, std::map & renumbering_map) { AKANTU_DEBUG_IN(); UInt * connectivity = list_nodes; UInt new_node_num = renumbering_map.size(); for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) { UInt & node = *connectivity; std::map::iterator it = renumbering_map.find(node); if(it == renumbering_map.end()) { UInt old_node = node; renumbering_map[old_node] = new_node_num; node = new_node_num; ++new_node_num; } else { node = it->second; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::purifyMesh(Mesh & mesh) { AKANTU_DEBUG_IN(); std::map renumbering_map; RemovedNodesEvent remove_nodes(mesh); Array & nodes_removed = remove_nodes.getList(); for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; Mesh::type_iterator it = mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined); for(; it != end; ++it) { ElementType type(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Array & connectivity_vect = mesh.getConnectivity(type, ghost_type); UInt nb_element(connectivity_vect.getSize()); UInt * connectivity = connectivity_vect.storage(); renumberNodesInConnectivity (connectivity, nb_element*nb_nodes_per_element, renumbering_map); } } Array & new_numbering = remove_nodes.getNewNumbering(); std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1)); std::map::iterator it = renumbering_map.begin(); std::map::iterator end = renumbering_map.end(); for (; it != end; ++it) { new_numbering(it->first) = it->second; } for (UInt i = 0; i < new_numbering.getSize(); ++i) { if(new_numbering(i) == UInt(-1)) nodes_removed.push_back(i); } mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } #if defined(AKANTU_COHESIVE_ELEMENT) /* -------------------------------------------------------------------------- */ void MeshUtils::insertCohesiveElements(Mesh & mesh, Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion, Array & doubled_nodes, Array & new_elements) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); if (updateFacetToDouble(mesh_facets, facet_insertion)) { if (spatial_dimension == 1) { doublePointFacet(mesh, mesh_facets, doubled_nodes); } else { doubleFacet(mesh, mesh_facets, spatial_dimension - 1, doubled_nodes, true); findSubfacetToDouble(mesh, mesh_facets); if (spatial_dimension == 2) { doubleSubfacet<2>(mesh, mesh_facets, doubled_nodes); } else if (spatial_dimension == 3) { doubleFacet(mesh, mesh_facets, 1, doubled_nodes, false); findSubfacetToDouble(mesh, mesh_facets); doubleSubfacet<3>(mesh, mesh_facets, doubled_nodes); } } updateCohesiveData(mesh, mesh_facets, new_elements); } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doubleNodes(Mesh & mesh, const Array & old_nodes, Array & doubled_nodes) { AKANTU_DEBUG_IN(); Array & position = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt old_nb_nodes = position.getSize(); UInt new_nb_nodes = old_nb_nodes + old_nodes.getSize(); UInt old_nb_doubled_nodes = doubled_nodes.getSize(); UInt new_nb_doubled_nodes = old_nb_doubled_nodes + old_nodes.getSize(); position.resize(new_nb_nodes); doubled_nodes.resize(new_nb_doubled_nodes); Array::iterator > position_begin = position.begin(spatial_dimension); for (UInt n = 0; n < old_nodes.getSize(); ++n) { UInt new_node = old_nb_nodes + n; /// store doubled nodes doubled_nodes(old_nb_doubled_nodes + n, 0) = old_nodes(n); doubled_nodes(old_nb_doubled_nodes + n, 1) = new_node; /// update position std::copy(position_begin + old_nodes(n), position_begin + old_nodes(n) + 1, position_begin + new_node); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::doubleFacet(Mesh & mesh, Mesh & mesh_facets, UInt facet_dimension, Array & doubled_nodes, bool facet_mode) { AKANTU_DEBUG_IN(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(facet_dimension, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(facet_dimension, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); if (nb_facet_to_double == 0) continue; ElementType type_subfacet = Mesh::getFacetType(type_facet); const UInt nb_subfacet_per_facet = Mesh::getNbFacetsPerElement(type_facet); GhostType gt_subfacet = _casper; Array > * f_to_subfacet = NULL; Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); UInt nb_nodes_per_facet = conn_facet.getNbComponent(); UInt old_nb_facet = conn_facet.getSize(); UInt new_nb_facet = old_nb_facet + nb_facet_to_double; conn_facet.resize(new_nb_facet); subfacet_to_facet.resize(new_nb_facet); UInt new_facet = old_nb_facet - 1; Element new_facet_el(type_facet, 0, gt_facet); Array::iterator > subfacet_to_facet_begin = subfacet_to_facet.begin(nb_subfacet_per_facet); Array::iterator > conn_facet_begin = conn_facet.begin(nb_nodes_per_facet); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); ++new_facet; /// adding a new facet by copying original one /// copy connectivity in new facet std::copy(conn_facet_begin + old_facet, conn_facet_begin + old_facet + 1, conn_facet_begin + new_facet); /// update subfacet_to_facet std::copy(subfacet_to_facet_begin + old_facet, subfacet_to_facet_begin + old_facet + 1, subfacet_to_facet_begin + new_facet); new_facet_el.element = new_facet; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet_per_facet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; if (gt_subfacet != subfacet.ghost_type) { gt_subfacet = subfacet.ghost_type; f_to_subfacet = & mesh_facets.getElementToSubelement(type_subfacet, subfacet.ghost_type); } /// update facet_to_subfacet array (*f_to_subfacet)(subfacet.element).push_back(new_facet_el); } } /// update facet_to_subfacet and _segment_3 facets if any if (!facet_mode) { updateSubfacetToFacet(mesh_facets, type_facet, gt_facet, true); updateFacetToSubfacet(mesh_facets, type_facet, gt_facet, true); updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } else updateQuadraticSegments(mesh, mesh_facets, type_facet, gt_facet, doubled_nodes); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool MeshUtils::updateFacetToDouble(Mesh & mesh_facets, const ElementTypeMapArray & facet_insertion) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); bool facets_to_double = false; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; const Array & f_insertion = facet_insertion(type_facet, gt_facet); Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); ElementType el_type = _not_defined; GhostType el_gt = _casper; UInt nb_facet_per_element = 0; Element old_facet_el(type_facet, 0, gt_facet); Array * facet_to_element = NULL; for (UInt f = 0; f < f_insertion.getSize(); ++f) { if (f_insertion(f) == false) continue; facets_to_double = true; if (element_to_facet(f)[1].type == _not_defined #if defined(AKANTU_COHESIVE_ELEMENT) || element_to_facet(f)[1].kind == _ek_cohesive #endif ) { AKANTU_DEBUG_WARNING("attempt to double a facet on the boundary"); continue; } f_to_double.push_back(f); UInt new_facet = mesh_facets.getNbElement(type_facet, gt_facet) + f_to_double.getSize() - 1; old_facet_el.element = f; /// update facet_to_element vector Element & elem_to_update = element_to_facet(f)[1]; UInt el = elem_to_update.element; if (elem_to_update.ghost_type != el_gt || elem_to_update.type != el_type) { el_type = elem_to_update.type; el_gt = elem_to_update.ghost_type; facet_to_element = & mesh_facets.getSubelementToElement(el_type, el_gt); nb_facet_per_element = facet_to_element->getNbComponent(); } Element * f_update = std::find(facet_to_element->storage() + el * nb_facet_per_element, facet_to_element->storage() + el * nb_facet_per_element + nb_facet_per_element, old_facet_el); AKANTU_DEBUG_ASSERT(facet_to_element->storage() + el * nb_facet_per_element != facet_to_element->storage() + el * nb_facet_per_element + nb_facet_per_element, "Facet not found"); f_update->element = new_facet; /// update elements connected to facet std::vector first_facet_list = element_to_facet(f); element_to_facet.push_back(first_facet_list); /// set new and original facets as boundary facets element_to_facet(new_facet)[0] = element_to_facet(new_facet)[1]; element_to_facet(f)[1] = ElementNull; element_to_facet(new_facet)[1] = ElementNull; } } } AKANTU_DEBUG_OUT(); return facets_to_double; } /* -------------------------------------------------------------------------- */ void MeshUtils::resetFacetToDouble(Mesh & mesh_facets) { AKANTU_DEBUG_IN(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh_facets.firstType(_all_dimensions, gt); Mesh::type_iterator end = mesh_facets.lastType(_all_dimensions, gt); for(; it != end; ++it) { ElementType type = *it; mesh_facets.getDataPointer("facet_to_double", type, gt, 1, false); mesh_facets.getDataPointer > ("facets_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer > ("elements_to_subfacet_double", type, gt, 1, false); mesh_facets.getDataPointer > ("subfacets_to_subsubfacet_double", type, gt, 1, false); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::findSubfacetToDouble(Mesh & mesh, Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); if (spatial_dimension == 1) return; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); if (nb_facet_to_double == 0) continue; ElementType type_subfacet = Mesh::getFacetType(type_facet); GhostType gt_subfacet = _casper; ElementType type_subsubfacet = Mesh::getFacetType(type_subfacet); GhostType gt_subsubfacet = _casper; Array * conn_subfacet = NULL; Array * sf_to_double = NULL; Array > * sf_to_subfacet_double = NULL; Array > * f_to_subfacet_double = NULL; Array > * el_to_subfacet_double = NULL; UInt nb_subfacet = Mesh::getNbFacetsPerElement(type_facet); UInt nb_subsubfacet; UInt nb_nodes_per_sf_el; if (subsubfacet_mode) { nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subsubfacet); nb_subsubfacet = Mesh::getNbFacetsPerElement(type_subfacet); } else nb_nodes_per_sf_el = Mesh::getNbNodesPerElement(type_subfacet); Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); Array * subsubfacet_to_subfacet = NULL; UInt old_nb_facet = subfacet_to_facet.getSize() - nb_facet_to_double; Element current_facet(type_facet, 0, gt_facet); std::vector element_list; std::vector facet_list; std::vector * subfacet_list; if (subsubfacet_mode) subfacet_list = new std::vector; /// map to filter subfacets Array< std::vector > * facet_to_subfacet = NULL; /// this is used to make sure that both new and old facets are /// checked UInt facets[2]; /// loop on every facet for (UInt f_index = 0; f_index < 2; ++f_index) { for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { facets[bool(f_index)] = f_to_double(facet); facets[!bool(f_index)] = old_nb_facet + facet; UInt old_facet = facets[0]; UInt new_facet = facets[1]; Element & starting_element = element_to_facet(new_facet)[0]; current_facet.element = old_facet; /// loop on every subfacet for (UInt sf = 0; sf < nb_subfacet; ++sf) { Element & subfacet = subfacet_to_facet(old_facet, sf); if (subfacet == ElementNull) continue; if (gt_subfacet != subfacet.ghost_type) { gt_subfacet = subfacet.ghost_type; if (subsubfacet_mode) { subsubfacet_to_subfacet = & mesh_facets.getSubelementToElement(type_subfacet, gt_subfacet); } else { conn_subfacet = & mesh_facets.getConnectivity(type_subfacet, gt_subfacet); sf_to_double = & mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); f_to_subfacet_double = & mesh_facets.getData >("facets_to_subfacet_double", type_subfacet, gt_subfacet); el_to_subfacet_double = & mesh_facets.getData >("elements_to_subfacet_double", type_subfacet, gt_subfacet); facet_to_subfacet = & mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet); } } if (subsubfacet_mode) { /// loop on every subsubfacet for (UInt ssf = 0; ssf < nb_subsubfacet; ++ssf) { Element & subsubfacet = (*subsubfacet_to_subfacet)(subfacet.element, ssf); if (subsubfacet == ElementNull) continue; if (gt_subsubfacet != subsubfacet.ghost_type) { gt_subsubfacet = subsubfacet.ghost_type; conn_subfacet = & mesh_facets.getConnectivity(type_subsubfacet, gt_subsubfacet); sf_to_double = & mesh_facets.getData("facet_to_double", type_subsubfacet, gt_subsubfacet); sf_to_subfacet_double = & mesh_facets.getData >("subfacets_to_subsubfacet_double", type_subsubfacet, gt_subsubfacet); f_to_subfacet_double = & mesh_facets.getData >("facets_to_subfacet_double", type_subsubfacet, gt_subsubfacet); el_to_subfacet_double = & mesh_facets.getData >("elements_to_subfacet_double", type_subsubfacet, gt_subsubfacet); facet_to_subfacet = & mesh_facets.getElementToSubelement(type_subsubfacet, gt_subsubfacet); } UInt global_ssf = subsubfacet.element; Vector subsubfacet_connectivity(conn_subfacet->storage() + global_ssf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subsubfacet is to be doubled if (findElementsAroundSubfacet(mesh, mesh_facets, starting_element, current_facet, subsubfacet_connectivity, element_list, facet_list, subfacet_list) == false && removeElementsInVector(*subfacet_list, (*facet_to_subfacet)(global_ssf)) == false) { sf_to_double->push_back(global_ssf); sf_to_subfacet_double->push_back(*subfacet_list); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } else { const UInt global_sf = subfacet.element; Vector subfacet_connectivity(conn_subfacet->storage() + global_sf * nb_nodes_per_sf_el, nb_nodes_per_sf_el); /// check if subfacet is to be doubled if (findElementsAroundSubfacet(mesh, mesh_facets, starting_element, current_facet, subfacet_connectivity, element_list, facet_list) == false && removeElementsInVector(facet_list, (*facet_to_subfacet)(global_sf)) == false) { sf_to_double->push_back(global_sf); f_to_subfacet_double->push_back(facet_list); el_to_subfacet_double->push_back(element_list); } } } } } if (subsubfacet_mode) delete subfacet_list; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) void MeshUtils::updateCohesiveData(Mesh & mesh, Mesh & mesh_facets, Array & new_elements) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); bool third_dimension = spatial_dimension == 3; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); if (nb_facet_to_double == 0) continue; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); Array * facet_to_coh_element = mesh_facets.getSubelementToElementPointer(type_cohesive, gt_facet); Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array * conn_cohesive = mesh.getConnectivityPointer(type_cohesive, gt_facet); UInt nb_nodes_per_facet = Mesh::getNbNodesPerElement(type_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); UInt old_nb_cohesive_elements = conn_cohesive->getSize(); UInt new_nb_cohesive_elements = conn_cohesive->getSize() + nb_facet_to_double; UInt old_nb_facet = element_to_facet.getSize() - nb_facet_to_double; facet_to_coh_element->resize(new_nb_cohesive_elements); conn_cohesive->resize(new_nb_cohesive_elements); UInt new_elements_old_size = new_elements.getSize(); new_elements.resize(new_elements_old_size + nb_facet_to_double); Element c_element(type_cohesive, 0, gt_facet, _ek_cohesive); Element f_element(type_facet, 0, gt_facet); UInt facets[2]; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { /// (in 3D cohesive elements connectivity is inverted) facets[third_dimension] = f_to_double(facet); facets[!third_dimension] = old_nb_facet + facet; UInt cohesive_element = old_nb_cohesive_elements + facet; /// store doubled facets f_element.element = facets[0]; (*facet_to_coh_element)(cohesive_element, 0) = f_element; f_element.element = facets[1]; (*facet_to_coh_element)(cohesive_element, 1) = f_element; /// modify cohesive elements' connectivity for (UInt n = 0; n < nb_nodes_per_facet; ++n) { (*conn_cohesive)(cohesive_element, n) = conn_facet(facets[0], n); (*conn_cohesive)(cohesive_element, n + nb_nodes_per_facet) = conn_facet(facets[1], n); } /// update element_to_facet vectors c_element.element = cohesive_element; element_to_facet(facets[0])[1] = c_element; element_to_facet(facets[1])[1] = c_element; /// add cohesive element to the element event list new_elements(new_elements_old_size + facet) = c_element; } } } AKANTU_DEBUG_OUT(); } #endif /* -------------------------------------------------------------------------- */ void MeshUtils::doublePointFacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); if (spatial_dimension != 1) return; Array & position = mesh.getNodes(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_facet = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); for(; it != end; ++it) { ElementType type_facet = *it; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); const Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); UInt old_nb_facet = element_to_facet.getSize() - nb_facet_to_double; UInt new_nb_facet = element_to_facet.getSize(); UInt old_nb_nodes = position.getSize(); UInt new_nb_nodes = old_nb_nodes + nb_facet_to_double; position.resize(new_nb_nodes); conn_facet.resize(new_nb_facet); UInt old_nb_doubled_nodes = doubled_nodes.getSize(); doubled_nodes.resize(old_nb_doubled_nodes + nb_facet_to_double); for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt new_facet = old_nb_facet + facet; ElementType type = element_to_facet(new_facet)[0].type; UInt el = element_to_facet(new_facet)[0].element; GhostType gt = element_to_facet(new_facet)[0].ghost_type; UInt old_node = conn_facet(old_facet); UInt new_node = old_nb_nodes + facet; /// update position position(new_node) = position(old_node); conn_facet(new_facet) = new_node; Array & conn_segment = mesh.getConnectivity(type, gt); UInt nb_nodes_per_segment = conn_segment.getNbComponent(); /// update facet connectivity UInt i; for (i = 0; conn_segment(el, i) != old_node && i <= nb_nodes_per_segment; ++i); conn_segment(el, i) = new_node; doubled_nodes(old_nb_doubled_nodes + facet, 0) = old_node; doubled_nodes(old_nb_doubled_nodes + facet, 1) = new_node; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::updateQuadraticSegments(Mesh & mesh, Mesh & mesh_facets, ElementType type_facet, GhostType gt_facet, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (type_facet != _segment_3) return; Array & f_to_double = mesh_facets.getData("facet_to_double", type_facet, gt_facet); UInt nb_facet_to_double = f_to_double.getSize(); UInt old_nb_facet = mesh_facets.getNbElement(type_facet, gt_facet) - nb_facet_to_double; Array & conn_facet = mesh_facets.getConnectivity(type_facet, gt_facet); Array< std::vector > & element_to_facet = mesh_facets.getElementToSubelement(type_facet, gt_facet); /// this ones matter only for segments in 3D Array< std::vector > * el_to_subfacet_double = NULL; Array< std::vector > * f_to_subfacet_double = NULL; if (third_dim_segments) { el_to_subfacet_double = & mesh_facets.getData >("elements_to_subfacet_double", type_facet, gt_facet); f_to_subfacet_double = & mesh_facets.getData >("facets_to_subfacet_double", type_facet, gt_facet); } Array middle_nodes; Array subfacets; for (UInt facet = 0; facet < nb_facet_to_double; ++facet) { UInt old_facet = f_to_double(facet); UInt node = conn_facet(old_facet, 2); if (!mesh.isPureGhostNode(node)) { middle_nodes.push_back(node); subfacets.push_back(facet); } } UInt old_nb_doubled_nodes = doubled_nodes.getSize(); doubleNodes(mesh, middle_nodes, doubled_nodes); for (UInt n = old_nb_doubled_nodes; n < doubled_nodes.getSize(); ++n) { UInt old_node = doubled_nodes(n, 0); UInt new_node = doubled_nodes(n, 1); UInt sf = subfacets(n - old_nb_doubled_nodes); UInt new_facet = old_nb_facet + sf; conn_facet(new_facet, 2) = new_node; if (third_dim_segments) { updateElementalConnectivity(mesh_facets, old_node, new_node, element_to_facet(new_facet)); updateElementalConnectivity(mesh, old_node, new_node, (*el_to_subfacet_double)(sf), &(*f_to_subfacet_double)(sf)); } else { updateElementalConnectivity(mesh, old_node, new_node, element_to_facet(new_facet)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateSubfacetToFacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.getSize(); /// update subfacet_to_facet vector ElementType type_facet = _not_defined; GhostType gt_facet = _casper; Array * subfacet_to_facet = NULL; UInt nb_subfacet_per_facet = 0; UInt old_nb_subfacet = mesh_facets.getNbElement(type_subfacet, gt_subfacet) - nb_subfacet_to_double; Array > * facet_list = NULL; if (facet_mode) facet_list = & mesh_facets.getData >("facets_to_subfacet_double", type_subfacet, gt_subfacet); else facet_list = & mesh_facets.getData >("subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); Element old_subfacet_el(type_subfacet, 0, gt_subfacet); Element new_subfacet_el(type_subfacet, 0, gt_subfacet); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { old_subfacet_el.element = sf_to_double(sf); new_subfacet_el.element = old_nb_subfacet + sf; for (UInt f = 0; f < (*facet_list)(sf).size(); ++f) { Element & facet = (*facet_list)(sf)[f]; if (facet.type != type_facet || facet.ghost_type != gt_facet) { type_facet = facet.type; gt_facet = facet.ghost_type; subfacet_to_facet = & mesh_facets.getSubelementToElement(type_facet, gt_facet); nb_subfacet_per_facet = subfacet_to_facet->getNbComponent(); } Element * sf_update = std::find(subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet, subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, old_subfacet_el); AKANTU_DEBUG_ASSERT(subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet != subfacet_to_facet->storage() + facet.element * nb_subfacet_per_facet + nb_subfacet_per_facet, "Subfacet not found"); *sf_update = new_subfacet_el; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::updateFacetToSubfacet(Mesh & mesh_facets, ElementType type_subfacet, GhostType gt_subfacet, bool facet_mode) { AKANTU_DEBUG_IN(); Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.getSize(); Array< std::vector > & facet_to_subfacet = mesh_facets.getElementToSubelement(type_subfacet, gt_subfacet); Array< std::vector > * facet_to_subfacet_double = NULL; if (facet_mode) { facet_to_subfacet_double = & mesh_facets.getData >("facets_to_subfacet_double", type_subfacet, gt_subfacet); } else { facet_to_subfacet_double = & mesh_facets.getData >("subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); } for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) facet_to_subfacet.push_back((*facet_to_subfacet_double)(sf)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshUtils::doubleSubfacet(Mesh & mesh, Mesh & mesh_facets, Array & doubled_nodes) { AKANTU_DEBUG_IN(); if (spatial_dimension == 1) return; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType gt_subfacet = *gt; Mesh::type_iterator it = mesh_facets.firstType(0, gt_subfacet); Mesh::type_iterator end = mesh_facets.lastType(0, gt_subfacet); for(; it != end; ++it) { ElementType type_subfacet = *it; Array & sf_to_double = mesh_facets.getData("facet_to_double", type_subfacet, gt_subfacet); UInt nb_subfacet_to_double = sf_to_double.getSize(); if (nb_subfacet_to_double == 0) continue; AKANTU_DEBUG_ASSERT(type_subfacet == _point_1, "Only _point_1 subfacet doubling is supported at the moment"); Array & conn_subfacet = mesh_facets.getConnectivity(type_subfacet, gt_subfacet); UInt old_nb_subfacet = conn_subfacet.getSize(); UInt new_nb_subfacet = old_nb_subfacet + nb_subfacet_to_double; conn_subfacet.resize(new_nb_subfacet); Array nodes_to_double; UInt old_nb_doubled_nodes = doubled_nodes.getSize(); /// double nodes for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_subfacet = sf_to_double(sf); nodes_to_double.push_back(conn_subfacet(old_subfacet)); } doubleNodes(mesh, nodes_to_double, doubled_nodes); /// add new nodes in connectivity for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt new_subfacet = old_nb_subfacet + sf; UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); conn_subfacet(new_subfacet) = new_node; } /// update facet and element connectivity Array > & f_to_subfacet_double = mesh_facets.getData >("facets_to_subfacet_double", type_subfacet, gt_subfacet); Array > & el_to_subfacet_double = mesh_facets.getData >("elements_to_subfacet_double", type_subfacet, gt_subfacet); Array > * sf_to_subfacet_double = NULL; if (spatial_dimension == 3) sf_to_subfacet_double = & mesh_facets.getData >("subfacets_to_subsubfacet_double", type_subfacet, gt_subfacet); for (UInt sf = 0; sf < nb_subfacet_to_double; ++sf) { UInt old_node = doubled_nodes(old_nb_doubled_nodes + sf, 0); UInt new_node = doubled_nodes(old_nb_doubled_nodes + sf, 1); updateElementalConnectivity(mesh, old_node, new_node, el_to_subfacet_double(sf), &f_to_subfacet_double(sf)); updateElementalConnectivity(mesh_facets, old_node, new_node, f_to_subfacet_double(sf)); if (spatial_dimension == 3) updateElementalConnectivity(mesh_facets, old_node, new_node, (*sf_to_subfacet_double)(sf)); } if (spatial_dimension == 2) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, true); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, true); } else if (spatial_dimension == 3) { updateSubfacetToFacet(mesh_facets, type_subfacet, gt_subfacet, false); updateFacetToSubfacet(mesh_facets, type_subfacet, gt_subfacet, false); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::flipFacets(Mesh & mesh_facets, const ElementTypeMapArray & global_connectivity, GhostType gt_facet) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh_facets.getSpatialDimension(); /// get global connectivity for local mesh ElementTypeMapArray global_connectivity_tmp; mesh_facets.initElementTypeMapArray(global_connectivity_tmp, 1, spatial_dimension - 1, gt_facet, true, _ek_regular, true); mesh_facets.getGlobalConnectivity(global_connectivity_tmp, spatial_dimension - 1, gt_facet); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, gt_facet); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, gt_facet); /// loop on every facet for(; it != end; ++it) { ElementType type_facet = *it; Array & connectivity = mesh_facets.getConnectivity(type_facet, gt_facet); const Array & g_connectivity = global_connectivity(type_facet, gt_facet); Array > & el_to_f = mesh_facets.getElementToSubelement(type_facet, gt_facet); Array & subfacet_to_facet = mesh_facets.getSubelementToElement(type_facet, gt_facet); UInt nb_subfacet_per_facet = subfacet_to_facet.getNbComponent(); UInt nb_nodes_per_facet = connectivity.getNbComponent(); UInt nb_facet = connectivity.getSize(); UInt nb_nodes_per_P1_facet = Mesh::getNbNodesPerElement(Mesh::getP1ElementType(type_facet)); Array & global_conn_tmp = global_connectivity_tmp(type_facet, gt_facet); Array::iterator > conn_it = connectivity.begin(nb_nodes_per_facet); Array::iterator > gconn_tmp_it = global_conn_tmp.begin(nb_nodes_per_facet); Array::const_iterator > conn_glob_it = g_connectivity.begin(nb_nodes_per_facet); Array::iterator > subf_to_f = subfacet_to_facet.begin(nb_subfacet_per_facet); UInt * conn_tmp_pointer = new UInt[nb_nodes_per_facet]; Vector conn_tmp(conn_tmp_pointer, nb_nodes_per_facet); Element el_tmp; Element * subf_tmp_pointer = new Element[nb_subfacet_per_facet]; Vector subf_tmp(subf_tmp_pointer, nb_subfacet_per_facet); for (UInt f = 0; f < nb_facet; ++f, ++conn_it, ++gconn_tmp_it, ++subf_to_f, ++conn_glob_it) { Vector & gconn_tmp = *gconn_tmp_it; const Vector & conn_glob = *conn_glob_it; Vector & conn_local = *conn_it; /// skip facet if connectivities are the same if (gconn_tmp == conn_glob) continue; /// re-arrange connectivity conn_tmp = conn_local; UInt * begin = conn_glob.storage(); UInt * end = conn_glob.storage() + nb_nodes_per_facet; for (UInt n = 0; n < nb_nodes_per_facet; ++n) { UInt * new_node = std::find(begin, end, gconn_tmp(n)); AKANTU_DEBUG_ASSERT(new_node != end, "Node not found"); UInt new_position = new_node - begin; conn_local(new_position) = conn_tmp(n); } /// if 3D, check if facets are just rotated if (spatial_dimension == 3) { /// find first node UInt * new_node = std::find(begin, end, gconn_tmp(0)); AKANTU_DEBUG_ASSERT(new_node != end, "Node not found"); UInt new_position = new_node - begin; UInt n = 1; /// count how many nodes in the received connectivity follow /// the same order of those in the local connectivity for (; n < nb_nodes_per_P1_facet && gconn_tmp(n) == conn_glob((new_position + n) % nb_nodes_per_P1_facet); ++n); /// skip the facet inversion if facet is just rotated if (n == nb_nodes_per_P1_facet) continue; } /// update data to invert facet el_tmp = el_to_f(f)[0]; el_to_f(f)[0] = el_to_f(f)[1]; el_to_f(f)[1] = el_tmp; subf_tmp = (*subf_to_f); (*subf_to_f)(0) = subf_tmp(1); (*subf_to_f)(1) = subf_tmp(0); } delete [] subf_tmp_pointer; delete [] conn_tmp_pointer; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::fillElementToSubElementsData(Mesh & mesh) { AKANTU_DEBUG_IN(); if(mesh.getNbElement(mesh.getSpatialDimension() - 1) == 0) { AKANTU_DEBUG_INFO("There are not facets, add them in the mesh file or call the buildFacet method."); return; } UInt spatial_dimension = mesh.getSpatialDimension(); ElementTypeMapArray barycenters("barycenter_tmp", mesh.getID()); mesh.initElementTypeMapArray(barycenters, spatial_dimension, _all_dimensions); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator it = mesh.firstType(_all_dimensions, *gt); Mesh::type_iterator end = mesh.lastType(_all_dimensions, *gt); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, *gt); Array & barycenters_arr = barycenters(*it, *gt); barycenters_arr.resize(nb_element); Array::vector_iterator bary = barycenters_arr.begin(spatial_dimension); Array::vector_iterator bary_end = barycenters_arr.end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { mesh.getBarycenter(el, *it, bary->storage(), *gt); } } } for(Int sp(spatial_dimension); sp >= 1; --sp) { if(mesh.getNbElement(sp) == 0) continue; for (ghost_type_t::iterator git = ghost_type_t::begin(); git != ghost_type_t::end(); ++git) { Mesh::type_iterator tit = mesh.firstType(sp, *git); Mesh::type_iterator tend = mesh.lastType(sp, *git); for (;tit != tend; ++tit) { mesh.getSubelementToElementPointer(*tit, *git)->resize(mesh.getNbElement(*tit, *git)); mesh.getSubelementToElement(*tit, *git).clear(); } tit = mesh.firstType(sp - 1, *git); tend = mesh.lastType(sp - 1, *git); for (;tit != tend; ++tit) { mesh.getElementToSubelementPointer(*tit, *git)->resize(mesh.getNbElement(*tit, *git)); mesh.getElementToSubelement(*tit, *git).clear(); } } CSR nodes_to_elements; buildNode2Elements(mesh, nodes_to_elements, sp); Element facet_element; for (ghost_type_t::iterator git = ghost_type_t::begin(); git != ghost_type_t::end(); ++git) { Mesh::type_iterator tit = mesh.firstType(sp - 1, *git); Mesh::type_iterator tend = mesh.lastType(sp - 1, *git); facet_element.ghost_type = *git; for (;tit != tend; ++tit) { facet_element.type = *tit; Array< std::vector > & element_to_subelement = mesh.getElementToSubelement(*tit, *git); const Array & connectivity = mesh.getConnectivity(*tit, *git); Array::const_iterator< Vector > fit = connectivity.begin(mesh.getNbNodesPerElement(*tit)); Array::const_iterator< Vector > fend = connectivity.end(mesh.getNbNodesPerElement(*tit)); UInt fid = 0; for (;fit != fend; ++fit, ++fid) { const Vector & facet = *fit; facet_element.element = fid; std::map element_seen_counter; UInt nb_nodes_per_facet = mesh.getNbNodesPerElement(Mesh::getP1ElementType(*tit)); for (UInt n(0); n < nb_nodes_per_facet; ++n) { CSR::iterator eit = nodes_to_elements.begin(facet(n)); CSR::iterator eend = nodes_to_elements.end(facet(n)); for(;eit != eend; ++eit) { Element & elem = *eit; std::map::iterator cit = element_seen_counter.find(elem); if(cit != element_seen_counter.end()) { cit->second++; } else { element_seen_counter[elem] = 1; } } } std::vector connected_elements; std::map::iterator cit = element_seen_counter.begin(); std::map::iterator cend = element_seen_counter.end(); for(;cit != cend; ++cit) { if(cit->second == nb_nodes_per_facet) connected_elements.push_back(cit->first); } std::vector::iterator ceit = connected_elements.begin(); std::vector::iterator ceend = connected_elements.end(); for(;ceit != ceend; ++ceit) element_to_subelement(fid).push_back(*ceit); for (UInt ce = 0; ce < connected_elements.size(); ++ce) { Element & elem = connected_elements[ce]; Array & subelement_to_element = *(mesh.getSubelementToElementPointer(elem.type, elem.ghost_type)); UInt f(0); for(; f < mesh.getNbFacetsPerElement(elem.type) && subelement_to_element(elem.element, f) != ElementNull; ++f); AKANTU_DEBUG_ASSERT(f < mesh.getNbFacetsPerElement(elem.type), "The element " << elem << " seems to have too many facets!! (" << f << " < " << mesh.getNbFacetsPerElement(elem.type) << ")"); subelement_to_element(elem.element, f) = facet_element; } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MeshUtils::findElementsAroundSubfacet(const Mesh & mesh, const Mesh & mesh_facets, const Element & starting_element, const Element & end_facet, const Vector & subfacet_connectivity, std::vector & elem_list, std::vector & facet_list, std::vector * subfacet_list) { AKANTU_DEBUG_IN(); /// preallocated stuff before starting bool facet_matched = false; elem_list.resize(0); facet_list.resize(0); if (third_dim_points) subfacet_list->resize(0); elem_list.push_back(starting_element); const Array * facet_connectivity = NULL; const Array * sf_connectivity = NULL; const Array * facet_to_element = NULL; const Array * subfacet_to_facet = NULL; ElementType current_type = _not_defined; GhostType current_ghost_type = _casper; ElementType current_facet_type = _not_defined; GhostType current_facet_ghost_type = _casper; ElementType current_subfacet_type = _not_defined; GhostType current_subfacet_ghost_type = _casper; const Array< std::vector > * element_to_facet = NULL; const Element * opposing_el = NULL; Array elements_to_check; elements_to_check.resize(0); elements_to_check.push_back(starting_element); /// keep going until there are elements to check while (elements_to_check.getSize() != 0) { /// loop over each element to check for (UInt el = 0; el < elements_to_check.getSize(); ++el) { Element & current_el = elements_to_check(el); if (current_el.type != current_type || current_el.ghost_type != current_ghost_type) { current_type = current_el.type; current_ghost_type = current_el.ghost_type; facet_to_element = & mesh_facets.getSubelementToElement(current_type, current_ghost_type); } /// loop over each facet of the element for (UInt f = 0; f < facet_to_element->getNbComponent(); ++f) { const Element & current_facet = (*facet_to_element)(current_el.element, f); if (current_facet == ElementNull) continue; if (current_facet_type != current_facet.type || current_facet_ghost_type != current_facet.ghost_type) { current_facet_type = current_facet.type; current_facet_ghost_type = current_facet.ghost_type; element_to_facet = & mesh_facets.getElementToSubelement(current_facet_type, current_facet_ghost_type); facet_connectivity = & mesh_facets.getConnectivity(current_facet_type, current_facet_ghost_type); if (third_dim_points) subfacet_to_facet = & mesh_facets.getSubelementToElement(current_facet_type, current_facet_ghost_type); } /// check if end facet is reached if (current_facet == end_facet) facet_matched = true; /// add this facet if not already passed if (std::find(facet_list.begin(), facet_list.end(), current_facet) == facet_list.end() && hasElement(*facet_connectivity, current_facet, subfacet_connectivity)) { facet_list.push_back(current_facet); if (third_dim_points) { /// check subfacets for (UInt sf = 0; sf < subfacet_to_facet->getNbComponent(); ++sf) { const Element & current_subfacet = (*subfacet_to_facet)(current_facet.element, sf); if (current_subfacet == ElementNull) continue; if (current_subfacet_type != current_subfacet.type || current_subfacet_ghost_type != current_subfacet.ghost_type) { current_subfacet_type = current_subfacet.type; current_subfacet_ghost_type = current_subfacet.ghost_type; sf_connectivity = & mesh_facets.getConnectivity(current_subfacet_type, current_subfacet_ghost_type); } if (std::find(subfacet_list->begin(), subfacet_list->end(), current_subfacet) == subfacet_list->end() && hasElement(*sf_connectivity, current_subfacet, subfacet_connectivity)) subfacet_list->push_back(current_subfacet); } } } else continue; /// consider opposing element if ( (*element_to_facet)(current_facet.element)[0] == current_el) opposing_el = & (*element_to_facet)(current_facet.element)[1]; else opposing_el = & (*element_to_facet)(current_facet.element)[0]; /// skip null elements since they are on a boundary if (*opposing_el == ElementNull) continue; /// skip this element if already added if ( std::find(elem_list.begin(), elem_list.end(), *opposing_el) != elem_list.end() ) continue; /// only regular elements have to be checked if (opposing_el->kind == _ek_regular) elements_to_check.push_back(*opposing_el); elem_list.push_back(*opposing_el); #ifndef AKANTU_NDEBUG const Array & conn_elem = mesh.getConnectivity(opposing_el->type, opposing_el->ghost_type); AKANTU_DEBUG_ASSERT(hasElement(conn_elem, *opposing_el, subfacet_connectivity), "Subfacet doesn't belong to this element"); #endif } /// erased checked element from the list elements_to_check.erase(el); } } AKANTU_DEBUG_OUT(); return facet_matched; } __END_AKANTU__ // LocalWords: ElementType diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index 3e4800bae..c9bf451cc 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,1569 +1,1652 @@ /** * @file material.cc * * @author Aurelia Isabel Cuba Ramos * @author Marco Vocialta * @author Nicolas Richart * @author Daniel Pino Muñoz * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Sep 16 2014 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #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(id, model.getMemoryID()), Parsable(_st_material, id), is_init(false), finite_deformation(false), name(""), model(&model), spatial_dimension(this->model->getSpatialDimension()), element_filter("element_filter", id, this->memory_id), stress("stress", *this), eigenstrain("eigenstrain", *this), gradu("grad_u", *this), green_strain("green_strain",*this), piola_kirchhoff_2("piola_kirchhoff_2", *this), // potential_energy_vector(false), potential_energy("potential_energy", *this), is_non_local(false), use_previous_stress(false), use_previous_gradu(false), interpolation_inverse_coordinates("interpolation inverse coordinates", *this), interpolation_points_matrices("interpolation points matrices", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the material model.getMesh().initElementTypeMapArray(element_filter, 1, spatial_dimension, false, _ek_regular); registerParam("rho" , rho , 0. , _pat_parsable | _pat_modifiable, "Density"); registerParam("name" , name , std::string(), _pat_parsable | _pat_readable); registerParam("finite_deformation" , finite_deformation , false , _pat_parsable | _pat_readable, "Is finite deformation"); registerParam("inelastic_deformation", inelastic_deformation, false , _pat_internal, "Is inelastic deformation"); /// allocate gradu stress for local elements eigenstrain.initialize(spatial_dimension * spatial_dimension); gradu.initialize(spatial_dimension * spatial_dimension); stress.initialize(spatial_dimension * spatial_dimension); model.registerEventHandler(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); if(finite_deformation) { this->piola_kirchhoff_2.initialize(spatial_dimension * spatial_dimension); if(use_previous_stress) this->piola_kirchhoff_2.initializeHistory(); this->green_strain.initialize(spatial_dimension * spatial_dimension); } if(use_previous_stress) this->stress.initializeHistory(); if(use_previous_gradu) this->gradu.initializeHistory(); for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::savePreviousState() { AKANTU_DEBUG_IN(); for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { if(it->second->hasHistory()) it->second->saveCurrentValues(); } 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(GhostType ghost_type) { AKANTU_DEBUG_IN(); computeAllStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); if(!finite_deformation){ Array & residual = const_cast &>(model->getResidual()); Mesh & mesh = model->getFEEngine().getMesh(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { const Array & shapes_derivatives = model->getFEEngine().getShapesDerivatives(*it, ghost_type); Array & 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->getFEEngine().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$ Array * sigma_dphi_dx = new Array(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Array * shapesd_filtered = new Array(0, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, *it, ghost_type, elem_filter); Array & stress_vect = stress(*it, ghost_type); Array::matrix_iterator sigma = stress_vect.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator B = shapesd_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_sigma_it = sigma_dphi_dx->begin(spatial_dimension, nb_nodes_per_element); for (UInt q = 0; q < nb_element*nb_quadrature_points; ++q, ++sigma, ++B, ++Bt_sigma_it) Bt_sigma_it->mul(*sigma, *B); 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$ */ Array * int_sigma_dphi_dx = new Array(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); model->getFEEngine().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, *it, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble model->getFEEngine().assembleArray(*int_sigma_dphi_dx, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, -1); delete int_sigma_dphi_dx; } } else{ switch (spatial_dimension){ case 1: this->assembleResidual<1>(ghost_type); break; case 2: this->assembleResidual<2>(ghost_type); break; case 3: this->assembleResidual<3>(ghost_type); break; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the gradu * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEEngine().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEEngine().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); Array & gradu_vect = gradu(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEEngine().gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect, spatial_dimension, *it, ghost_type, elem_filter); gradu_vect -= eigenstrain(*it, ghost_type); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllCauchyStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(finite_deformation,"The Cauchy stress can only be computed if you are working in finite deformation."); //resizeInternalArray(stress); Mesh::type_iterator it = model->getFEEngine().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEEngine().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) switch (spatial_dimension){ case 1: this->computeCauchyStress<1>(*it, ghost_type); break; case 2: this->computeCauchyStress<2>(*it, ghost_type); break; case 3: this->computeCauchyStress<3>(*it, ghost_type); break; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeCauchyStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Array::matrix_iterator gradu_it = this->gradu(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator gradu_end = this->gradu(el_type, ghost_type).end(dim, dim); Array::matrix_iterator piola_it = this->piola_kirchhoff_2(el_type, ghost_type).begin(dim, dim); Array::matrix_iterator stress_it = this->stress(el_type, ghost_type).begin(dim, dim); Matrix F_tensor(dim, dim); for (; gradu_it != gradu_end; ++gradu_it, ++piola_it, ++stress_it) { Matrix & grad_u = *gradu_it; Matrix & piola = *piola_it; Matrix & sigma = *stress_it; gradUToF (grad_u, F_tensor); this->computeCauchyStressOnQuad(F_tensor, piola, sigma); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & displacement = model->getDisplacement(); //resizeInternalArray(gradu); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEEngine().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEEngine().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); Array & gradu_vect = gradu(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEEngine().gradientOnQuadraturePoints(displacement, gradu_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(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { if(finite_deformation){ switch (spatial_dimension) { case 1: { assembleStiffnessMatrixNL < 1 > (*it, ghost_type); assembleStiffnessMatrixL2 < 1 > (*it, ghost_type); break; } case 2: { assembleStiffnessMatrixNL < 2 > (*it, ghost_type); assembleStiffnessMatrixL2 < 2 > (*it, ghost_type); break; } case 3: { assembleStiffnessMatrixNL < 3 > (*it, ghost_type); assembleStiffnessMatrixL2 < 3 > (*it, ghost_type); break; } } } else { switch(spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(*it, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(*it, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(*it, ghost_type); break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); const Array & shapes_derivatives = model->getFEEngine().getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); model->getFEEngine().gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array(nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapesd_filtered = new Array(0, dim * nb_nodes_per_element, "filtered shapesd"); FEEngine::filterElementalData(model->getFEEngine().getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element, dim*nb_nodes_per_element); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::matrix_iterator D_end = tangent_stiffness_matrix->end (tangent_size, tangent_size); for(; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it) { Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix( *shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_D.mul(B, D); Bt_D_B.mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); model->getFEEngine().integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast (model->getStiffnessMatrix()); const Array & shapes_derivatives = model->getFEEngine().getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); //Array & gradu_vect = delta_gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(type, ghost_type); //gradu_vect.resize(nb_quadrature_points * nb_element); // model->getFEEngine().gradientOnQuadraturePoints(model->getIncrement(), gradu_vect, // dim, type, ghost_type, &elem_filter); Array * shapes_derivatives_filtered = new Array (nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_s_b_size = dim * nb_nodes_per_element; Array * bt_s_b = new Array (nb_element * nb_quadrature_points, bt_s_b_size * bt_s_b_size, "B^t*D*B"); UInt piola_matrix_size = getCauchyStressMatrixSize(dim); Matrix B(piola_matrix_size, bt_s_b_size); Matrix Bt_S(bt_s_b_size, piola_matrix_size); Matrix S(piola_matrix_size, piola_matrix_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_S_B_it = bt_s_b->begin(bt_s_b_size, bt_s_b_size); Array::matrix_iterator Bt_S_B_end = bt_s_b->end(bt_s_b_size, bt_s_b_size); Array::matrix_iterator piola_it = piola_kirchhoff_2(type, ghost_type).begin(dim, dim); for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) { Matrix & Bt_S_B = *Bt_S_B_it; Matrix & Piola_kirchhoff_matrix = *piola_it; setCauchyStressMatrix< dim >(Piola_kirchhoff_matrix, S); VoigtHelper::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_S.mul < true, false > (B, S); Bt_S_B.mul < false, false > (Bt_S, B); } delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array (nb_element, bt_s_b_size * bt_s_b_size, "K_e"); model->getFEEngine().integrate(*bt_s_b, *K_e, bt_s_b_size * bt_s_b_size, type, ghost_type, elem_filter); delete bt_s_b; model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast (model->getStiffnessMatrix()); const Array & shapes_derivatives = model->getFEEngine().getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); model->getFEEngine().gradientOnQuadraturePoints(model->getDisplacement(), gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Array * tangent_stiffness_matrix = new Array (nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Array * shapes_derivatives_filtered = new Array (nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Array::const_matrix_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Array * bt_d_b = new Array (nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); Matrix B(tangent_size, dim * nb_nodes_per_element); Matrix B2(tangent_size, dim * nb_nodes_per_element); Matrix Bt_D(dim * nb_nodes_per_element, tangent_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); Array::matrix_iterator Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element, dim * nb_nodes_per_element); Array::matrix_iterator grad_u_it = gradu_vect.begin(dim, dim); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Array::matrix_iterator D_end = tangent_stiffness_matrix->end(tangent_size, tangent_size); for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) { Matrix & grad_u = *grad_u_it; Matrix & D = *D_it; Matrix & Bt_D_B = *Bt_D_B_it; //transferBMatrixToBL1 (*shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element); B += B2; Bt_D.mul < true, false > (B, D); Bt_D_B.mul < false, false > (Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * K_e = new Array (nb_element, bt_d_b_size * bt_d_b_size, "K_e"); model->getFEEngine().integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, elem_filter); delete bt_d_b; model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleResidual(GhostType ghost_type){ AKANTU_DEBUG_IN(); Array & residual = const_cast &> (model->getResidual()); Mesh & mesh = model->getFEEngine().getMesh(); Mesh::type_iterator it = element_filter.firstType(dim, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(dim, ghost_type); for (; it != last_type; ++it) { const Array & shapes_derivatives = model->getFEEngine().getShapesDerivatives(*it, ghost_type); Array & elem_filter = element_filter(*it, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(*it, ghost_type); Array * shapesd_filtered = new Array(0, size_of_shapes_derivatives, "filtered shapesd"); FEEngine::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered, *it, ghost_type, elem_filter); Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); //Set stress vectors UInt stress_size = getTangentStiffnessVoigtSize(dim); //Set matrices B and BNL* UInt bt_s_size = dim * nb_nodes_per_element; Array * bt_s = new Array (nb_element * nb_quadrature_points, bt_s_size, "B^t*S"); Array::matrix_iterator grad_u_it = this->gradu(*it, ghost_type).begin(dim, dim); Array::matrix_iterator grad_u_end = this->gradu(*it, ghost_type).end(dim, dim); Array::matrix_iterator stress_it = this->piola_kirchhoff_2(*it, ghost_type).begin(dim, dim); shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator bt_s_it = bt_s->begin(bt_s_size, 1); Matrix S_vect(stress_size, 1); Matrix B_tensor(stress_size, bt_s_size); Matrix B2_tensor(stress_size, bt_s_size); for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) { Matrix & grad_u = *grad_u_it; Matrix & r_it = *bt_s_it; Matrix & S_it = *stress_it; SetCauchyStressArray (S_it, S_vect); VoigtHelper::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element); VoigtHelper::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element); B_tensor += B2_tensor; r_it.mul < true, false > (B_tensor, S_vect); } delete shapesd_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Array * r_e = new Array (nb_element, bt_s_size, "r_e"); model->getFEEngine().integrate(*bt_s, *r_e, bt_s_size, *it, ghost_type, elem_filter); delete bt_s; model->getFEEngine().assembleArray(*r_e, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, -1); delete r_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { switch(spatial_dimension) { case 1: { computeAllStressesFromTangentModuli<1>(*it, ghost_type); break; } case 2: { computeAllStressesFromTangentModuli<2>(*it, ghost_type); break; } case 3: { computeAllStressesFromTangentModuli<3>(*it, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Array & shapes_derivatives = model->getFEEngine().getShapesDerivatives(type, ghost_type); Array & elem_filter = element_filter(type, ghost_type); Array & gradu_vect = gradu(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(type, ghost_type); gradu_vect.resize(nb_quadrature_points * nb_element); Array & disp = model->getDisplacement(); model->getFEEngine().gradientOnQuadraturePoints(disp, gradu_vect, dim, type, ghost_type, elem_filter); UInt tangent_moduli_size = getTangentStiffnessVoigtSize(dim); Array * tangent_moduli_tensors = new Array(nb_element*nb_quadrature_points, tangent_moduli_size * tangent_moduli_size, "tangent_moduli_tensors"); tangent_moduli_tensors->clear(); computeTangentModuli(type, *tangent_moduli_tensors, ghost_type); Array * shapesd_filtered = new Array(0, dim* nb_nodes_per_element, "filtered shapesd"); FEEngine::filterElementalData(model->getFEEngine().getMesh(), shapes_derivatives, *shapesd_filtered, type, ghost_type, elem_filter); Array filtered_u(nb_element, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(model->getFEEngine().getMesh(), disp, filtered_u, type, ghost_type, elem_filter); /// compute @f$\mathbf{D} \mathbf{B} \mathbf{u}@f$ Array::matrix_iterator shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element); Array::matrix_iterator D_it = tangent_moduli_tensors->begin(tangent_moduli_size, tangent_moduli_size); Array::matrix_iterator sigma_it = stress(type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::vector_iterator u_it = filtered_u.begin(spatial_dimension * nb_nodes_per_element); Matrix B(tangent_moduli_size, spatial_dimension * nb_nodes_per_element); Vector Bu(tangent_moduli_size); Vector DBu(tangent_moduli_size); for (UInt e = 0; e < nb_element; ++e, ++u_it) { for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it, ++shapes_derivatives_filtered_it, ++sigma_it) { Vector & u = *u_it; Matrix & sigma = *sigma_it; Matrix & D = *D_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bu.mul(B, u); DBu.mul(D, Bu); // Voigt notation to full symmetric tensor for (UInt i = 0; i < dim; ++i) sigma(i, i) = DBu(i); if(dim == 2) { sigma(0,1) = sigma(1,0) = DBu(2); } else if(dim == 3) { sigma(1,2) = sigma(2,1) = DBu(3); sigma(0,2) = sigma(2,0) = DBu(4); sigma(0,1) = sigma(1,0) = DBu(5); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElements() { AKANTU_DEBUG_IN(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension); for(; it != last_type; ++it) { computePotentialEnergy(*it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(!potential_energy.exists(el_type, ghost_type)) { UInt nb_element = element_filter(el_type, ghost_type).getSize(); UInt nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(el_type, _not_ghost); potential_energy.alloc(nb_element * nb_quadrature_points, 1, el_type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElements(); /// integrate the potential energy for each type of elements Mesh::type_iterator it = element_filter.firstType(spatial_dimension); Mesh::type_iterator last_type = element_filter.lastType(spatial_dimension); for(; it != last_type; ++it) { epot += model->getFEEngine().integrate(potential_energy(*it, _not_ghost), *it, _not_ghost, element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy(ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real epot = 0.; Vector epot_on_quad_points(model->getFEEngine().getNbQuadraturePoints(type)); computePotentialEnergyByElement(type, index, epot_on_quad_points); epot = model->getFEEngine().integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if(type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if(energy_id == "potential") return getPotentialEnergy(type, index); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::computeQuadraturePointsCoordinates(ElementTypeMapArray & quadrature_points_coordinates, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); - const Mesh & mesh = this->model->getFEEngine().getMesh(); + const Mesh & mesh = this->model->getMesh(); Array nodes_coordinates(mesh.getNodes(), true); nodes_coordinates += this->model->getDisplacement(); Mesh::type_iterator it = this->element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { const Array & elem_filter = this->element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_tot_quad = this->model->getFEEngine().getNbQuadraturePoints(*it, ghost_type) * nb_element; Array & quads = quadrature_points_coordinates(*it, ghost_type); quads.resize(nb_tot_quad); this->model->getFEEngine().interpolateOnQuadraturePoints(nodes_coordinates, quads, spatial_dimension, *it, ghost_type, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation(const ElementTypeMapArray & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEEngine().getMesh(); ElementTypeMapArray quadrature_points_coordinates("quadrature_points_coordinates_for_interpolation", getID()); mesh.initElementTypeMapArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; computeQuadraturePointsCoordinates(quadrature_points_coordinates, ghost_type); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last = element_filter.lastType(spatial_dimension, ghost_type); for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; const Array & interp_points_coord = interpolation_points_coordinates(type, ghost_type); UInt nb_interpolation_points_per_elem = interp_points_coord.getSize() / nb_element; AKANTU_DEBUG_ASSERT(interp_points_coord.getSize() % nb_element == 0, "Number of interpolation points is wrong"); #define AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD(type) \ initElementalFieldInterpolation(quadrature_points_coordinates(type, ghost_type), \ interp_points_coord, \ nb_interpolation_points_per_elem, \ ghost_type) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD); #undef AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::initElementalFieldInterpolation(const Array & quad_coordinates, const Array & interpolation_points_coordinates, const UInt nb_interpolation_points_per_elem, const GhostType ghost_type) { AKANTU_DEBUG_IN(); - UInt size_inverse_coords = - getSizeElementalFieldInterpolationCoodinates(ghost_type); - Array & elem_fil = element_filter(type, ghost_type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEEngine().getNbQuadraturePoints(type, ghost_type); if(!interpolation_inverse_coordinates.exists(type, ghost_type)) interpolation_inverse_coordinates.alloc(nb_element, - size_inverse_coords*size_inverse_coords, + nb_quad_per_element*nb_quad_per_element, type, ghost_type); if(!interpolation_points_matrices.exists(type, ghost_type)) interpolation_points_matrices.alloc(nb_element, - nb_interpolation_points_per_elem * size_inverse_coords, + nb_interpolation_points_per_elem * nb_quad_per_element, type, ghost_type); Array & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type); Array & interp_points_mat = interpolation_points_matrices(type, ghost_type); - Matrix quad_coord_matrix(size_inverse_coords, size_inverse_coords); + Matrix quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); Array::const_matrix_iterator quad_coords_it = quad_coordinates.begin_reinterpret(spatial_dimension, nb_quad_per_element, nb_element); Array::const_matrix_iterator points_coords_begin = interpolation_points_coordinates.begin_reinterpret(spatial_dimension, nb_interpolation_points_per_elem, interpolation_points_coordinates.getSize() / nb_interpolation_points_per_elem); Array::matrix_iterator inv_quad_coord_it = - interp_inv_coord.begin(size_inverse_coords, size_inverse_coords); + interp_inv_coord.begin(nb_quad_per_element, nb_quad_per_element); Array::matrix_iterator inv_points_mat_it = - interp_points_mat.begin(nb_interpolation_points_per_elem, size_inverse_coords); + interp_points_mat.begin(nb_interpolation_points_per_elem, nb_quad_per_element); /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el, ++inv_quad_coord_it, ++inv_points_mat_it, ++quad_coords_it) { /// matrix containing the quadrature points coordinates const Matrix & quad_coords = *quad_coords_it; /// matrix to store the matrix inversion result 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); /// matrix containing the interpolation points coordinates const Matrix & points_coords = points_coords_begin[elem_fil(el)]; /// matrix to store the interpolation points coordinates /// compatible with these functions Matrix & inv_points_coord_matrix = *inv_points_mat_it; /// insert the quad coordinates in a matrix compatible with the interpolation buildElementalFieldInterpolationCoodinates(points_coords, inv_points_coord_matrix); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type) { - AKANTU_DEBUG_IN(); + interpolateElementalField(stress, result, ghost_type); +} - const Mesh & mesh = model->getFEEngine().getMesh(); +/* -------------------------------------------------------------------------- */ +void Material::interpolateElementalField(const ElementTypeMapArray & field, + ElementTypeMapArray & result, + const GhostType ghost_type) { + AKANTU_DEBUG_IN(); Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last = element_filter.lastType(spatial_dimension, ghost_type); for (; it != last; ++it) { ElementType type = *it; - UInt nb_element = mesh.getNbElement(type, ghost_type); - if (nb_element == 0) continue; - Array & res = result(type, ghost_type); + Array & elem_fil = element_filter(type, ghost_type); + UInt nb_element = elem_fil.getSize(); + UInt nb_quad_per_element = model->getFEEngine().getNbQuadraturePoints(type, ghost_type); + + const Array & field_vec = field(type, ghost_type); + Array & result_vec = result(type, ghost_type); + + Matrix coefficients(nb_quad_per_element, field_vec.getNbComponent()); + + const Array & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type); + const Array & interp_points_coord = interpolation_points_matrices(type, ghost_type); + + UInt nb_interpolation_points_per_elem + = interp_points_coord.getNbComponent() / nb_quad_per_element; + + Array::const_matrix_iterator field_it + = field_vec.begin_reinterpret(field_vec.getNbComponent(), + nb_quad_per_element, + nb_element); + + Array::const_matrix_iterator interpolation_points_coordinates_it = + interp_points_coord.begin(nb_interpolation_points_per_elem, nb_quad_per_element); + + Array::matrix_iterator result_begin + = result_vec.begin_reinterpret(field_vec.getNbComponent(), + nb_interpolation_points_per_elem, + result_vec.getSize() / nb_interpolation_points_per_elem); -#define INTERPOLATE_ELEMENTAL_FIELD(type) \ - interpolateElementalField(stress(type, ghost_type), \ - res, \ - ghost_type) + Array::const_matrix_iterator inv_quad_coord_it = + interp_inv_coord.begin(nb_quad_per_element, nb_quad_per_element); - AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTERPOLATE_ELEMENTAL_FIELD); -#undef INTERPOLATE_ELEMENTAL_FIELD + /// loop over the elements of the current material and element type + for (UInt el = 0; el < nb_element; + ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { + /** + * matrix containing the inversion of the quadrature points' + * coordinates + */ + const 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); + + /// matrix containing the points' coordinates + const Matrix & coord = *interpolation_points_coordinates_it; + + /// multiply the coordinates matrix by the coefficients matrix and store the result + result_begin[elem_fil(el)].mul(coefficients, coord); + } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void Material::interpolateElementalField(const Array & field, - Array & result, +void Material::interpolateStressOnFacets(ElementTypeMapArray & result, const GhostType ghost_type) { + interpolateElementalFieldOnFacets(stress, result, ghost_type); +} + +/* -------------------------------------------------------------------------- */ +void Material::interpolateElementalFieldOnFacets(const ElementTypeMapArray & field, + ElementTypeMapArray & result, + const GhostType ghost_type) { AKANTU_DEBUG_IN(); - Array & elem_fil = element_filter(type, ghost_type); - UInt nb_element = elem_fil.getSize(); - UInt nb_quad_per_element = model->getFEEngine().getNbQuadraturePoints(type, ghost_type); - UInt size_inverse_coords = getSizeElementalFieldInterpolationCoodinates(ghost_type); + UInt sp2 = spatial_dimension * spatial_dimension; - Matrix coefficients(nb_quad_per_element, field.getNbComponent()); + const Mesh & mesh = this->model->getMesh(); + const Mesh & mesh_facets = mesh.getMeshFacets(); - const Array & interp_inv_coord = interpolation_inverse_coordinates(type, - ghost_type); - const Array & interp_points_coord = interpolation_points_matrices(type, - ghost_type); + Mesh::type_iterator it = element_filter.firstType(spatial_dimension, ghost_type); + Mesh::type_iterator last = element_filter.lastType(spatial_dimension, ghost_type); + for (; it != last; ++it) { + ElementType type = *it; - UInt nb_interpolation_points_per_elem = interp_points_coord.getNbComponent() / size_inverse_coords; + Array & elem_fil = element_filter(type, ghost_type); + UInt nb_element = elem_fil.getSize(); + UInt nb_quad_per_element = model->getFEEngine().getNbQuadraturePoints(type, ghost_type); - Array::const_matrix_iterator field_it - = field.begin_reinterpret(field.getNbComponent(), - nb_quad_per_element, - nb_element); + const Array & field_vec = field(type, ghost_type); - Array::const_matrix_iterator interpolation_points_coordinates_it = - interp_points_coord.begin(nb_interpolation_points_per_elem, size_inverse_coords); + Matrix coefficients(nb_quad_per_element, field_vec.getNbComponent()); - Array::matrix_iterator result_begin - = result.begin_reinterpret(field.getNbComponent(), - nb_interpolation_points_per_elem, - result.getSize() / nb_interpolation_points_per_elem); + const Array & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type); + const Array & interp_points_coord = interpolation_points_matrices(type, ghost_type); - Array::const_matrix_iterator inv_quad_coord_it = - interp_inv_coord.begin(size_inverse_coords, size_inverse_coords); + UInt nb_interpolation_points_per_elem + = interp_points_coord.getNbComponent() / nb_quad_per_element; - /// loop over the elements of the current material and element type - for (UInt el = 0; el < nb_element; - ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { - /** - * matrix containing the inversion of the quadrature points' - * coordinates - */ - const 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); - - /// matrix containing the points' coordinates - const Matrix & coord = *interpolation_points_coordinates_it; - - /// multiply the coordinates matrix by the coefficients matrix and store the result - result_begin[elem_fil(el)].mul(coefficients, coord); - } + Array::const_matrix_iterator field_it + = field_vec.begin_reinterpret(field_vec.getNbComponent(), + nb_quad_per_element, + nb_element); + + Array::const_matrix_iterator interpolation_points_coordinates_it = + interp_points_coord.begin(nb_interpolation_points_per_elem, nb_quad_per_element); + + Array::const_matrix_iterator inv_quad_coord_it = + interp_inv_coord.begin(nb_quad_per_element, nb_quad_per_element); + + Matrix result_tmp(sp2, nb_interpolation_points_per_elem); + const Array & facet_to_element = + mesh_facets.getSubelementToElement(type, ghost_type); + ElementType type_facet = Mesh::getFacetType(type); + UInt nb_facet_per_elem = facet_to_element.getNbComponent(); + UInt nb_quad_per_facet = nb_interpolation_points_per_elem / nb_facet_per_elem; + Element element_for_comparison(type, 0, ghost_type); + + /// loop over the elements of the current material and element type + for (UInt el = 0; el < nb_element; + ++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { + /** + * matrix containing the inversion of the quadrature points' + * coordinates + */ + const 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); + + /// matrix containing the points' coordinates + const Matrix & coord = *interpolation_points_coordinates_it; + + /// multiply the coordinates matrix by the coefficients matrix and store the result + result_tmp.mul(coefficients, coord); + + UInt global_el = elem_fil(el); + element_for_comparison.element = global_el; + for (UInt f = 0; f < nb_facet_per_elem; ++f) { + Element facet_elem = facet_to_element(global_el, f); + UInt global_facet = facet_elem.element; + const Array< std::vector > & element_to_facet + = mesh_facets.getElementToSubelement(type_facet, facet_elem.ghost_type); + Array & result_vec = result(type_facet, facet_elem.ghost_type); + + bool is_second_element = element_to_facet(global_facet)[0] != element_for_comparison; + + for (UInt q = 0; q < nb_quad_per_facet; ++q) { + Vector result_local(result_vec.storage() + + (global_facet * nb_quad_per_facet + q) * result_vec.getNbComponent() + + is_second_element * sp2, + sp2); + + result_local = result_tmp(f * nb_quad_per_facet + q); + } + } + } + } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const Array & Material::getArray(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 << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch(debug::Exception & e) { AKANTU_EXCEPTION("The material " << name << "(" < & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << getID() << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getArray(fvect_id); } catch(debug::Exception & e) { AKANTU_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]"); } } /* -------------------------------------------------------------------------- */ const InternalField & Material::getInternal(const ID & int_id) const { std::map *>::const_iterator it = internal_vectors_real.find(getID() + ":" + int_id); if(it == internal_vectors_real.end()) { AKANTU_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ InternalField & Material::getInternal(const ID & int_id) { std::map *>::iterator it = internal_vectors_real.find(getID() + ":" + int_id); if(it == internal_vectors_real.end()) { AKANTU_EXCEPTION("The material " << name << "(" << getID() << ") does not contain an internal " << int_id << " (" << (getID() + ":" + int_id) << ")"); } return *it->second; } /* -------------------------------------------------------------------------- */ void Material::addElements(const Array & elements_to_add) { AKANTU_DEBUG_IN(); UInt mat_id = model->getInternalIndexFromID(getID()); Array::const_iterator el_begin = elements_to_add.begin(); Array::const_iterator el_end = elements_to_add.end(); for(;el_begin != el_end; ++el_begin) { const Element & element = *el_begin; Array & mat_indexes = model->getMaterialByElement (element.type, element.ghost_type); Array & mat_loc_num = model->getMaterialLocalNumbering(element.type, element.ghost_type); UInt index = this->addElement(element.type, element.element, element.ghost_type); mat_indexes(element.element) = mat_id; mat_loc_num(element.element) = index; } this->resizeInternals(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::removeElements(const Array & elements_to_remove) { AKANTU_DEBUG_IN(); Array::const_iterator el_begin = elements_to_remove.begin(); Array::const_iterator el_end = elements_to_remove.end(); if(el_begin==el_end) return; ElementTypeMapArray material_local_new_numbering("remove mat filter elem", getID()); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; element.ghost_type = ghost_type; ElementTypeMapArray::type_iterator it = element_filter.firstType(_all_dimensions, ghost_type, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, ghost_type, _ek_not_defined); for(; it != end; ++it) { ElementType type = *it; element.type = type; Array & elem_filter = this->element_filter(type, ghost_type); Array & mat_loc_num = this->model->getMaterialLocalNumbering(type, ghost_type); if(!material_local_new_numbering.exists(type, ghost_type)) material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, ghost_type); Array & mat_renumbering = material_local_new_numbering(type, ghost_type); UInt nb_element = elem_filter.getSize(); element.kind=(*el_begin).kind; Array elem_filter_tmp; UInt new_id = 0; for (UInt el = 0; el < nb_element; ++el) { element.element = elem_filter(el); if(std::find(el_begin, el_end, element) == el_end) { elem_filter_tmp.push_back(element.element); mat_renumbering(el) = new_id; mat_loc_num(element.element) = new_id; ++new_id; } else { mat_renumbering(el) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.getSize()); elem_filter.copy(elem_filter_tmp); } } for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resizeInternals() { AKANTU_DEBUG_IN(); for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::onElementsAdded(__attribute__((unused)) const Array & element_list, __attribute__((unused)) const NewElementsEvent & event) { this->resizeInternals(); } /* -------------------------------------------------------------------------- */ void Material::onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { UInt my_num = model->getInternalIndexFromID(getID()); ElementTypeMapArray material_local_new_numbering("remove mat filter elem", getID()); Array::const_iterator el_begin = element_list.begin(); Array::const_iterator el_end = element_list.end(); for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType gt = *g; ElementTypeMapArray::type_iterator it = new_numbering.firstType(_all_dimensions, gt, _ek_not_defined); ElementTypeMapArray::type_iterator end = new_numbering.lastType (_all_dimensions, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if(element_filter.exists(type, gt)){ Array & elem_filter = element_filter(type, gt); Array & mat_indexes = this->model->getMaterialByElement (*it, gt); Array & mat_loc_num = this->model->getMaterialLocalNumbering(*it, gt); UInt nb_element = this->model->getMesh().getNbElement(type, gt); // all materials will resize of the same size... mat_indexes.resize(nb_element); mat_loc_num.resize(nb_element); if(!material_local_new_numbering.exists(type, gt)) material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, gt); Array & mat_renumbering = material_local_new_numbering(type, gt); const Array & renumbering = new_numbering(type, gt); Array elem_filter_tmp; UInt ni = 0; Element el; el.type = type; el.ghost_type = gt; for (UInt i = 0; i < elem_filter.getSize(); ++i) { el.element = elem_filter(i); if(std::find(el_begin, el_end, el) == el_end) { UInt new_el = renumbering(el.element); AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A not removed element as been badly renumbered"); elem_filter_tmp.push_back(new_el); mat_renumbering(i) = ni; mat_indexes(new_el) = my_num; mat_loc_num(new_el) = ni; ++ni; } else { mat_renumbering(i) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.getSize()); elem_filter.copy(elem_filter); } } } for (std::map *>::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); for (std::map *>::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->removeQuadraturePoints(material_local_new_numbering); } /* -------------------------------------------------------------------------- */ void Material::onBeginningSolveStep(const AnalysisMethod & method) { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onEndSolveStep(const AnalysisMethod & method) { ElementTypeMapArray::type_iterator it = this->element_filter.firstType(_all_dimensions, _not_ghost, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined); for(; it != end; ++it) { this->updateEnergies(*it, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDamageIteration() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void Material::onDamageUpdate() { ElementTypeMapArray::type_iterator it = this->element_filter.firstType(_all_dimensions, _not_ghost, _ek_not_defined); ElementTypeMapArray::type_iterator end = element_filter.lastType(_all_dimensions, _not_ghost, _ek_not_defined); for(; it != end; ++it) { if(!this->potential_energy.exists(*it, _not_ghost)) { UInt nb_element = this->element_filter(*it, _not_ghost).getSize(); UInt nb_quadrature_points = this->model->getFEEngine().getNbQuadraturePoints(*it, _not_ghost); this->potential_energy.alloc(nb_element * nb_quadrature_points, 1, *it, _not_ghost); } this->updateEnergiesAfterDamage(*it, _not_ghost); } } /* -------------------------------------------------------------------------- */ void Material::onDump(){ if(this->isFiniteDeformation()) this->computeAllCauchyStresses(_not_ghost); } /* -------------------------------------------------------------------------- */ void Material::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); std::string type = getID().substr(getID().find_last_of(":") + 1); stream << space << "Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Material::flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type, ElementKind element_kind){ typedef ElementTypeMapArray::type_iterator iterator; iterator tit = this->element_filter.firstType(this->spatial_dimension, ghost_type, element_kind); iterator end = this->element_filter.lastType(this->spatial_dimension, ghost_type, element_kind); for (; tit != end; ++tit) { ElementType type = *tit; try { __attribute__((unused)) const Array & src_vect = this->getArray(field_id,type,ghost_type); } catch(debug::Exception & e) { continue; } const Array & src_vect = this->getArray(field_id,type,ghost_type); const Array & filter = this->element_filter(type,ghost_type); // total number of elements for a given type UInt nb_element = this->model->mesh.getNbElement(type,ghost_type); // number of filtered elements UInt nb_element_src = filter.getSize(); // number of quadrature points per elem UInt nb_quad_per_elem = 0; // number of data per quadrature point UInt nb_data_per_quad = src_vect.getNbComponent(); if (!internal_flat.exists(type,ghost_type)) { internal_flat.alloc(nb_element*nb_quad_per_elem,nb_data_per_quad,type,ghost_type); } if (nb_element_src == 0) continue; nb_quad_per_elem = (src_vect.getSize()/nb_element_src); // number of data per element UInt nb_data = nb_quad_per_elem * src_vect.getNbComponent(); Array & dst_vect = internal_flat(type,ghost_type); dst_vect.resize(nb_element*nb_quad_per_elem); Array::const_scalar_iterator it = filter.begin(); Array::const_scalar_iterator end = filter.end(); Array::const_vector_iterator it_src = src_vect.begin_reinterpret(nb_data,nb_element_src); Array::vector_iterator it_dst = dst_vect.begin_reinterpret(nb_data,nb_element); for (; it != end ; ++it,++it_src) { it_dst[*it] = *it_src; } } }; /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 46186b53a..b53e92ae8 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,588 +1,596 @@ /** * @file material.hh * * @author Marco Vocialta * @author Nicolas Richart * @author Daniel Pino Muñoz * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Sep 16 2014 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_voigthelper.hh" #include "parser.hh" #include "parsable.hh" #include "data_accessor.hh" #include "internal_field.hh" #include "random_internal_field.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; } __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, * Array & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : public Memory, public DataAccessor, public Parsable, public MeshEventHandler, protected SolidMechanicsModelEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const ID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Function that materials can/should reimplement */ /* ------------------------------------------------------------------------ */ protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the tangent stiffness matrix virtual void computeTangentModuli(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the potential energy for an element virtual void computePotentialEnergyByElement(__attribute__((unused)) ElementType type, __attribute__((unused)) UInt index, __attribute__((unused)) Vector & epot_on_quad_points) { AKANTU_DEBUG_TO_IMPLEMENT(); } virtual void updateEnergies(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { } virtual void updateEnergiesAfterDamage(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} /// 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) { } /// function called to update the internal parameters when the modifiable /// parameters are modified virtual void updateInternalParameters() {} public: /// compute the p-wave speed in the material virtual Real getPushWaveSpeed(const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the s-wave speed in the material virtual Real getShearWaveSpeed(const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// get a material celerity to compute the stable time step (default: is the push wave speed) virtual Real getCelerity(const Element & element) const { return getPushWaveSpeed(element); } /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void registerInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } template void unregisterInternal(__attribute__((unused)) InternalField & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(GhostType ghost_type = _not_ghost); /// assemble the residual for this material virtual void assembleResidual(GhostType ghost_type); /// Operations before and after solveStep in implicit virtual void beforeSolveStep() {} virtual void afterSolveStep() {} /// save the stress in the previous_stress if needed virtual void savePreviousState(); /// compute the stresses for this material virtual void computeAllStresses(GhostType ghost_type = _not_ghost); virtual void computeAllNonLocalStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) {}; virtual void computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost); virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); /// add many elements at once void addElements(const Array & elements_to_add); /// remove many element at once void removeElements(const Array & elements_to_remove); /// 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 */ void interpolateStress(ElementTypeMapArray & result, const GhostType ghost_type = _not_ghost); + /** + * interpolate stress on given positions for each element by means + * of a geometrical interpolation on quadrature points and store the + * results per facet + */ + void interpolateStressOnFacets(ElementTypeMapArray & result, + const GhostType ghost_type = _not_ghost); + /** * function to initialize the elemental field interpolation * function by inverting the quadrature points' coordinates */ void initElementalFieldInterpolation(const ElementTypeMapArray & interpolation_points_coordinates); /* ------------------------------------------------------------------------ */ /* Common part */ /* ------------------------------------------------------------------------ */ protected: /// assemble the residual template void assembleResidual(GhostType ghost_type); /// Computation of Cauchy stress tensor in the case of finite deformation template void computeCauchyStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost); template inline void computeCauchyStressOnQuad(const Matrix & F, const Matrix & S, Matrix & cauchy, const Real & C33 = 1.0 ) const; template void computeAllStressesFromTangentModuli(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// assembling in finite deformation template void assembleStiffnessMatrixNL(const ElementType & type, GhostType ghost_type); template void assembleStiffnessMatrixL2(const ElementType & type, GhostType ghost_type); /// write the stress tensor in the Voigt notation. template inline void SetCauchyStressArray(const Matrix & S_t, Matrix & Stress_vect); inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// Size of the Stress matrix for the case of finite deformation see: Bathe et al, IJNME, Vol 9, 353-386, 1975 inline UInt getCauchyStressMatrixSize(UInt spatial_dimension) const; /// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386, 1975 template inline void setCauchyStressMatrix(const Matrix & S_t, Matrix & Stress_matrix); /// compute the potential energy by element void computePotentialEnergyByElements(); /// resize the intenals arrays void resizeInternals(); public: /// compute the coordinates of the quadrature points void computeQuadraturePointsCoordinates(ElementTypeMapArray & quadrature_points_coordinates, const GhostType & ghost_type) const; protected: /// interpolate an elemental field on given points for each element - template - void interpolateElementalField(const Array & field, - Array & result, + void interpolateElementalField(const ElementTypeMapArray & field, + ElementTypeMapArray & result, const GhostType ghost_type); + /// interpolate an elemental field and store the results per facet + void interpolateElementalFieldOnFacets(const ElementTypeMapArray & field, + ElementTypeMapArray & result, + const GhostType ghost_type); + /// template function to initialize the elemental field interpolation template void initElementalFieldInterpolation(const Array & quad_coordinates, const Array & interpolation_points_coordinates, const UInt nb_interpolation_points_per_elem, const GhostType ghost_type); /// build the coordinate matrix for the interpolation on elemental field template inline void buildElementalFieldInterpolationCoodinates(const Matrix & coordinates, Matrix & coordMatrix); /// build interpolation coordinates for basic linear elements inline void buildElementalFieldInterpolationCoodinatesLinear(const Matrix & coordinates, Matrix & coordMatrix); /// build interpolation coordinates for basic quadratic elements inline void buildElementalFieldInterpolationCoodinatesQuadratic(const Matrix & coordinates, Matrix & coordMatrix); - /// get the size of the coordiante matrix used in the interpolation - template - inline UInt getSizeElementalFieldInterpolationCoodinates(GhostType ghost_type = _not_ghost); - public: /* ------------------------------------------------------------------------ */ /* Conversion functions */ /* ------------------------------------------------------------------------ */ template inline void gradUToF (const Matrix & grad_u, Matrix & F) const; inline void rightCauchy(const Matrix & F, Matrix & C) const; inline void leftCauchy (const Matrix & F, Matrix & B) const; template inline void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) const; template inline void gradUToGreenStrain(const Matrix & grad_u, Matrix & epsilon) const; protected: /// converts global element to local element inline Element convertToLocalElement(const Element & global_element) const; /// converts local element to global element inline Element convertToGlobalElement(const Element & local_element) const; /// converts global quadrature point to local quadrature point inline QuadraturePoint convertToLocalPoint(const QuadraturePoint & global_point) const; /// converts local quadrature point to global quadrature point inline QuadraturePoint convertToGlobalPoint(const QuadraturePoint & local_point) const; /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; virtual inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; virtual inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); template inline void packElementDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); /* ------------------------------------------------------------------------ */ /* MeshEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ virtual void onElementsAdded(const Array & element_list, const NewElementsEvent & event); virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* SolidMechanicsModelEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: virtual void onBeginningSolveStep(const AnalysisMethod & method); virtual void onEndSolveStep(const AnalysisMethod & method); virtual void onDamageIteration(); virtual void onDamageUpdate(); virtual void onDump(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, Memory::getID(), const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); /// return the potential energy for the subset of elements contained by the material Real getPotentialEnergy(); /// return the potential energy for the provided element Real getPotentialEnergy(ElementType & type, UInt index); /// return the energy (identified by id) for the subset of elements contained by the material virtual Real getEnergy(std::string energy_id); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray &); AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray &); AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); bool isNonLocal() const { return is_non_local; } const Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; Array & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); const InternalField & getInternal(const ID & id) const; InternalField & getInternal(const ID & id); inline bool isInternal(const ID & id, const ElementKind & element_kind) const; inline ElementTypeMap getInternalDataPerElem(const ID & id, const ElementKind & element_kind) const; bool isFiniteDeformation() const { return finite_deformation; } bool isInelasticDeformation() const { return inelastic_deformation; } template inline void setParam(const ID & param, T value); template inline const T & getParam(const ID & param) const; void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, const GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined); protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// boolean to know if the material has been initialized bool is_init; std::map *> internal_vectors_real; std::map *> internal_vectors_uint; protected: /// Finite deformation bool finite_deformation; /// Finite deformation bool inelastic_deformation; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// spatial dimension UInt spatial_dimension; /// list of element handled by the material ElementTypeMapArray element_filter; /// stresses arrays ordered by element types InternalField stress; /// eigenstrain arrays ordered by element types InternalField eigenstrain; /// grad_u arrays ordered by element types InternalField gradu; /// Green Lagrange strain (Finite deformation) InternalField green_strain; /// Second Piola-Kirchhoff stress tensor arrays ordered by element types (Finite deformation) InternalField piola_kirchhoff_2; /// potential energy by element InternalField potential_energy; /// tell if using in non local mode or not bool is_non_local; /// tell if the material need the previous stress state bool use_previous_stress; /// tell if the material need the previous strain state bool use_previous_gradu; /// elemental field interpolation coordinates InternalField interpolation_inverse_coordinates; /// elemental field interpolation points InternalField interpolation_points_matrices; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ Array::matrix_iterator gradu_it = \ this->gradu(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ Array::matrix_iterator gradu_end = \ this->gradu(el_type, ghost_type).end(this->spatial_dimension, \ this->spatial_dimension); \ \ this->stress(el_type, \ ghost_type).resize(this->gradu(el_type, \ ghost_type).getSize()); \ \ Array::iterator< Matrix > stress_it = \ this->stress(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ \ if(this->isFiniteDeformation()){ \ this->piola_kirchhoff_2(el_type, \ ghost_type).resize(this->gradu(el_type, \ ghost_type).getSize()); \ stress_it = \ this->piola_kirchhoff_2(el_type, \ ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ } \ \ for(;gradu_it != gradu_end; ++gradu_it, ++stress_it) { \ Matrix & __attribute__((unused)) grad_u = *gradu_it; \ Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Array::matrix_iterator gradu_it = \ this->gradu(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ Array::matrix_iterator gradu_end = \ this->gradu(el_type, ghost_type).end(this->spatial_dimension, \ this->spatial_dimension); \ Array::matrix_iterator sigma_it = \ this->stress(el_type, ghost_type).begin(this->spatial_dimension, \ this->spatial_dimension); \ \ tangent_mat.resize(this->gradu(el_type, ghost_type).getSize()); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(this->spatial_dimension); \ Array::matrix_iterator tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;gradu_it != gradu_end; ++gradu_it, ++sigma_it, ++tangent_it) { \ Matrix & __attribute__((unused)) grad_u = *gradu_it; \ Matrix & __attribute__((unused)) sigma_tensor = *sigma_it; \ Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } \ /* -------------------------------------------------------------------------- */ #define INSTANSIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #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 f7e3fb1e4..9819833f6 100644 --- a/src/model/solid_mechanics/material_inline_impl.cc +++ b/src/model/solid_mechanics/material_inline_impl.cc @@ -1,482 +1,476 @@ /** * @file material_inline_impl.cc * * @author Marco Vocialta * @author Nicolas Richart * @author Daniel Pino Muñoz * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Sep 16 2014 * * @brief Implementation of the inline functions of the class material * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ __END_AKANTU__ #include "solid_mechanics_model.hh" #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const ElementType & type, UInt element, const GhostType & ghost_type) { Array & el_filter = element_filter(type, ghost_type); el_filter.push_back(element); return el_filter.getSize()-1; } /* -------------------------------------------------------------------------- */ inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const { return (dim * (dim - 1) / 2 + dim); } /* -------------------------------------------------------------------------- */ inline UInt Material::getCauchyStressMatrixSize(UInt dim) const { return (dim * dim); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToF(const Matrix & grad_u, Matrix & F) const { AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim*dim, "The dimension of the tensor F should be greater or equal to the dimension of the tensor grad_u."); F.eye(); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) F(i, j) += grad_u(i, j); } /* -------------------------------------------------------------------------- */ template inline void Material::computeCauchyStressOnQuad(const Matrix & F, const Matrix & piola, Matrix & sigma, const Real & C33 ) const { Real J = F.det() * sqrt(C33); Matrix F_S(dim, dim); F_S.mul(F, piola); Real constant = J ? 1./J : 0; sigma.mul(F_S, F, constant); } /* -------------------------------------------------------------------------- */ inline void Material::rightCauchy(const Matrix & F, Matrix & C) const { C.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::leftCauchy(const Matrix & F, Matrix & B) const { B.mul(F, F); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) const { for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) epsilon(i, j) = 0.5*(grad_u(i, j) + grad_u(j, i)); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToGreenStrain(const Matrix & grad_u, Matrix & epsilon) const { epsilon.mul(grad_u, grad_u, .5); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) epsilon(i, j) += 0.5 * (grad_u(i, j) + grad_u(j, i)); } /* ---------------------------------------------------------------------------*/ template inline void Material::SetCauchyStressArray(const Matrix & S_t, Matrix & Stress_vect) { AKANTU_DEBUG_IN(); Stress_vect.clear(); //UInt cauchy_matrix_size = getCauchyStressArraySize(dim); //see Finite ekement formulations for large deformation dynamic analysis, Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau /* * 1d: [ s11 ]' * 2d: [ s11 s22 s12 ]' * 3d: [ s11 s22 s33 s23 s13 s12 ] */ for (UInt i = 0; i < dim; ++i)//diagonal terms Stress_vect(i, 0) = S_t(i, i); for (UInt i = 1; i < dim; ++i)// term s12 in 2D and terms s23 s13 in 3D Stress_vect(dim+i-1, 0) = S_t(dim-i-1, dim-1); for (UInt i = 2; i < dim; ++i)//term s13 in 3D Stress_vect(dim+i, 0) = S_t(0, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void Material::setCauchyStressMatrix(const Matrix & S_t, Matrix & Stress_matrix) { AKANTU_DEBUG_IN(); Stress_matrix.clear(); /// see Finite ekement formulations for large deformation dynamic analysis, /// Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau for (UInt i = 0; i < dim; ++i) { for (UInt m = 0; m < dim; ++m) { for (UInt n = 0; n < dim; ++n) { Stress_matrix(i * dim + m, i * dim + n) = S_t(m, n); } } } //other terms from the diagonal /*for (UInt i = 0; i < 3 - dim; ++i) { Stress_matrix(dim * dim + i, dim * dim + i) = S_t(dim + i, dim + i); }*/ AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline Element Material::convertToLocalElement(const Element & global_element) const { UInt ge = global_element.element; #ifndef AKANTU_NDEBUG UInt model_mat_index = this->model->getMaterialByElement(global_element.type, global_element.ghost_type)(ge); UInt mat_index = this->model->getMaterialIndex(this->name); AKANTU_DEBUG_ASSERT(model_mat_index == mat_index, "Conversion of a global element in a local element for the wrong material " << this->name << std::endl); #endif UInt le = this->model->getMaterialLocalNumbering(global_element.type, global_element.ghost_type)(ge); Element tmp_quad(global_element.type, le, global_element.ghost_type); return tmp_quad; } /* -------------------------------------------------------------------------- */ inline Element Material::convertToGlobalElement(const Element & local_element) const { UInt le = local_element.element; UInt ge = this->element_filter(local_element.type, local_element.ghost_type)(le); Element tmp_quad(local_element.type, ge, local_element.ghost_type); return tmp_quad; } /* -------------------------------------------------------------------------- */ inline QuadraturePoint Material::convertToLocalPoint(const QuadraturePoint & global_point) const { const FEEngine & fem = this->model->getFEEngine(); UInt nb_quad = fem.getNbQuadraturePoints(global_point.type); Element el = this->convertToLocalElement(static_cast(global_point)); QuadraturePoint tmp_quad(el, global_point.num_point, nb_quad); return tmp_quad; } /* -------------------------------------------------------------------------- */ inline QuadraturePoint Material::convertToGlobalPoint(const QuadraturePoint & local_point) const { const FEEngine & fem = this->model->getFEEngine(); UInt nb_quad = fem.getNbQuadraturePoints(local_point.type); Element el = this->convertToGlobalElement(static_cast(local_point)); QuadraturePoint tmp_quad(el, local_point.num_point, nb_quad); return tmp_quad; } /* -------------------------------------------------------------------------- */ template inline void Material::buildElementalFieldInterpolationCoodinates(__attribute__((unused)) const Matrix & coordinates, __attribute__((unused)) Matrix & coordMatrix) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ inline void Material::buildElementalFieldInterpolationCoodinatesLinear(const Matrix & coordinates, Matrix & coordMatrix) { for (UInt i = 0; i < coordinates.cols(); ++i) coordMatrix(i, 0) = 1; } /* -------------------------------------------------------------------------- */ inline void Material::buildElementalFieldInterpolationCoodinatesQuadratic(const Matrix & coordinates, Matrix & coordMatrix) { UInt nb_quadrature_points = coordMatrix.cols(); for (UInt i = 0; i < coordinates.cols(); ++i) { coordMatrix(i, 0) = 1; for (UInt j = 1; j < nb_quadrature_points; ++j) coordMatrix(i, j) = coordinates(j-1, i); } } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_segment_2>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesLinear(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_segment_3>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesQuadratic(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_3>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesLinear(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_6>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesQuadratic(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_tetrahedron_4>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesLinear(coordinates, coordMatrix); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_tetrahedron_10>(const Matrix & coordinates, Matrix & coordMatrix) { buildElementalFieldInterpolationCoodinatesQuadratic(coordinates, coordMatrix); } /** * @todo Write a more efficient interpolation for quadrangles by * dropping unnecessary quadrature points * */ /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_4>(const Matrix & coordinates, Matrix & coordMatrix) { for (UInt i = 0; i < coordinates.cols(); ++i) { Real x = coordinates(0, i); Real y = coordinates(1, i); coordMatrix(i, 0) = 1; coordMatrix(i, 1) = x; coordMatrix(i, 2) = y; coordMatrix(i, 3) = x * y; } } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_8>(const Matrix & coordinates, Matrix & coordMatrix) { for (UInt i = 0; i < coordinates.cols(); ++i) { UInt j = 0; Real x = coordinates(0, i); Real y = coordinates(1, i); for (UInt e = 0; e <= 2; ++e) { for (UInt n = 0; n <= 2; ++n) { coordMatrix(i, j) = std::pow(x, e) * std::pow(y, n); ++j; } } } } -/* -------------------------------------------------------------------------- */ -template -inline UInt Material::getSizeElementalFieldInterpolationCoodinates(GhostType ghost_type) { - return model->getFEEngine().getNbQuadraturePoints(type, ghost_type); -} - /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataForElements(const Array & elements, SynchronizationTag tag) const { if(tag == _gst_smm_stress) { return (this->isFiniteDeformation() ? 3 : 1) * spatial_dimension * spatial_dimension * sizeof(Real) * this->getModel().getNbQuadraturePoints(elements); } return 0; } /* -------------------------------------------------------------------------- */ inline void Material::packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const { if(tag == _gst_smm_stress) { if(this->isFiniteDeformation()) { packElementDataHelper(piola_kirchhoff_2, buffer, elements); packElementDataHelper(gradu, buffer, elements); } packElementDataHelper(stress, buffer, elements); } } /* -------------------------------------------------------------------------- */ inline void Material::unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) { if(tag == _gst_smm_stress) { if(this->isFiniteDeformation()) { unpackElementDataHelper(piola_kirchhoff_2, buffer, elements); unpackElementDataHelper(gradu, buffer, elements); } unpackElementDataHelper(stress, buffer, elements); } } /* -------------------------------------------------------------------------- */ template inline const T & Material::getParam(const ID & param) const { try { return get(param); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID()); } } /* -------------------------------------------------------------------------- */ template inline void Material::setParam(const ID & param, T value) { try { set(param, value); } catch(...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID()); } updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template inline void Material::packElementDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) const { DataAccessor::packElementalDataHelper(data_to_pack, buffer, elements, true, model->getFEEngine(fem_id)); } /* -------------------------------------------------------------------------- */ template inline void Material::unpackElementDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) { DataAccessor::unpackElementalDataHelper(data_to_unpack, buffer, elements, true, model->getFEEngine(fem_id)); } /* -------------------------------------------------------------------------- */ template<> inline void Material::registerInternal(InternalField & vect) { internal_vectors_real[vect.getID()] = &vect; } template<> inline void Material::registerInternal(InternalField & vect) { internal_vectors_uint[vect.getID()] = &vect; } /* -------------------------------------------------------------------------- */ template<> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_real.erase(vect.getID()); } template<> inline void Material::unregisterInternal(InternalField & vect) { internal_vectors_uint.erase(vect.getID()); } /* -------------------------------------------------------------------------- */ inline bool Material::isInternal(const ID & id, const ElementKind & element_kind) const { std::map *>::const_iterator internal_array = internal_vectors_real.find(this->getID()+":"+id); if (internal_array == internal_vectors_real.end()) return false; if (internal_array->second->getElementKind() != element_kind) return false; return true; } /* -------------------------------------------------------------------------- */ inline ElementTypeMap Material::getInternalDataPerElem(const ID & id, const ElementKind & element_kind) const { std::map *>::const_iterator internal_array = internal_vectors_real.find(this->getID()+":"+id); if (internal_array == internal_vectors_real.end()) AKANTU_EXCEPTION("cannot find internal " << id); if (internal_array->second->getElementKind() != element_kind) AKANTU_EXCEPTION("cannot find internal " << id); InternalField & internal = *internal_array->second; InternalField::type_iterator it = internal.firstType(spatial_dimension, _not_ghost,element_kind); InternalField::type_iterator last_type = internal.lastType(spatial_dimension, _not_ghost,element_kind); ElementTypeMap res; for(; it != last_type; ++it) { UInt nb_quadrature_points = 0; if (element_kind == _ek_regular) nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(*it); #if defined(AKANTU_COHESIVE_ELEMENT) else if (element_kind == _ek_cohesive) nb_quadrature_points = model->getFEEngine("CohesiveFEEngine").getNbQuadraturePoints(*it); #endif res(*it) = internal.getNbComponent() * nb_quadrature_points; } return res; } diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc index 7d5c3ec56..6e59fa986 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc @@ -1,517 +1,467 @@ /** * @file material_cohesive_linear.cc * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Thu Aug 07 2014 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinear::MaterialCohesiveLinear(SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model,id), sigma_c_eff("sigma_c_eff", *this), delta_c_eff("delta_c_eff", *this), insertion_stress("insertion_stress", *this) { AKANTU_DEBUG_IN(); this->registerParam("beta" , beta , 0. , _pat_parsable | _pat_readable, "Beta parameter" ); this->registerParam("G_c" , G_c , 0. , _pat_parsable | _pat_readable, "Mode I fracture energy" ); this->registerParam("penalty", penalty, 0. , _pat_parsable | _pat_readable, "Penalty coefficient" ); this->registerParam("volume_s", volume_s, 0. , _pat_parsable | _pat_readable, "Reference volume for sigma_c scaling"); this->registerParam("m_s", m_s, 1. , _pat_parsable | _pat_readable, "Weibull exponent for sigma_c scaling"); this->registerParam("kappa" , kappa , 1. , _pat_parsable | _pat_readable, "Kappa parameter"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); /// compute scalars beta2_kappa2 = beta * beta/kappa/kappa; beta2_kappa = beta * beta/kappa; if (Math::are_float_equal(beta, 0)) beta2_inv = 0; else beta2_inv = 1./beta/beta; sigma_c_eff.initialize(1); delta_c_eff.initialize(1); insertion_stress.initialize(spatial_dimension); if (!Math::are_float_equal(delta_c, 0.)) delta_c_eff.setDefaultValue(delta_c); else delta_c_eff.setDefaultValue(2 * G_c / sigma_c); if (model->getIsExtrinsic()) scaleInsertionTraction(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::scaleInsertionTraction() { AKANTU_DEBUG_IN(); // do nothing if volume_s hasn't been specified by the user if (Math::are_float_equal(volume_s, 0.)) return; const Mesh & mesh_facets = model->getMeshFacets(); const FEEngine & fe_engine = model->getFEEngine(); const FEEngine & fe_engine_facet = model->getFEEngine("FacetsFEEngine"); // loop over facet type Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); Real base_sigma_c = sigma_c; for(;first != last; ++first) { ElementType type_facet = *first; const Array< std::vector > & facet_to_element = mesh_facets.getElementToSubelement(type_facet); UInt nb_facet = facet_to_element.getSize(); UInt nb_quad_per_facet = fe_engine_facet.getNbQuadraturePoints(type_facet); // iterator to modify sigma_c for all the quadrature points of a facet Array::vector_iterator sigma_c_iterator = sigma_c(type_facet).begin_reinterpret(nb_quad_per_facet, nb_facet); for (UInt f = 0; f < nb_facet; ++f, ++sigma_c_iterator) { const std::vector & element_list = facet_to_element(f); // compute bounding volume Real volume = 0; std::vector::const_iterator elem = element_list.begin(); std::vector::const_iterator elem_end = element_list.end(); for (; elem != elem_end; ++elem) { if (*elem == ElementNull) continue; // unit vector for integration in order to obtain the volume UInt nb_quadrature_points = fe_engine.getNbQuadraturePoints(elem->type); Vector unit_vector(nb_quadrature_points, 1); volume += fe_engine.integrate(unit_vector, elem->type, elem->element, elem->ghost_type); } // scale sigma_c *sigma_c_iterator -= base_sigma_c; *sigma_c_iterator *= std::pow(volume_s / volume, 1. / m_s); *sigma_c_iterator += base_sigma_c; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::checkInsertion() { AKANTU_DEBUG_IN(); const Mesh & mesh_facets = model->getMeshFacets(); CohesiveElementInserter & inserter = model->getElementInserter(); Real tolerance = Math::getTolerance(); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); for (; it != last; ++it) { ElementType type_facet = *it; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); const Array & facets_check = inserter.getCheckFacets(type_facet); Array & f_insertion = inserter.getInsertionFacets(type_facet); Array & f_filter = facet_filter(type_facet); Array & sig_c_eff = sigma_c_eff(type_cohesive); Array & del_c = delta_c_eff(type_cohesive); Array & ins_stress = insertion_stress(type_cohesive); Array & trac_old = tractions_old(type_cohesive); const Array & f_stress = model->getStressOnFacets(type_facet); const Array & sigma_lim = sigma_c(type_facet); - UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbQuadraturePoints(type_facet); + UInt nb_quad_facet = model->getFacetFEEngine().getNbQuadraturePoints(type_facet); UInt nb_facet = f_filter.getSize(); if (nb_facet == 0) continue; - UInt tot_nb_quad = nb_facet * nb_quad_facet; + Array::const_iterator sigma_lim_it = sigma_lim.begin(); - Array stress_check (tot_nb_quad); - Array normal_traction(tot_nb_quad, spatial_dimension); + Matrix stress_tmp(spatial_dimension, spatial_dimension); + Matrix normal_traction(spatial_dimension, nb_quad_facet); + Vector stress_check(nb_quad_facet); + UInt sp2 = spatial_dimension * spatial_dimension; + + const Array & tangents = model->getTangents(type_facet); + const Array & normals + = model->getFacetFEEngine().getNormalsOnQuadPoints(type_facet); + Array::const_vector_iterator normal_begin = normals.begin(spatial_dimension); + Array::const_vector_iterator tangent_begin = tangents.begin(tangents.getNbComponent()); + Array::const_matrix_iterator facet_stress_begin = + f_stress.begin(spatial_dimension, spatial_dimension * 2); + + std::list new_sigmas; + std::list< Vector > new_normal_traction; + std::list new_delta_c; + + // loop over each facet belonging to this material + for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it) { + UInt facet = f_filter(f); + // skip facets where check shouldn't be realized + if (!facets_check(facet)) continue; - computeStressNorms(f_stress, stress_check, normal_traction, type_facet); + // compute the effective norm on each quadrature point of the facet + for (UInt q = 0; q < nb_quad_facet; ++q) { + UInt current_quad = f * nb_quad_facet + q; + const Vector & normal = normal_begin[current_quad]; + const Vector & tangent = tangent_begin[current_quad]; + const Matrix & facet_stress_it = facet_stress_begin[current_quad]; - Array::iterator > stress_check_it = - stress_check.begin_reinterpret(nb_quad_facet, nb_facet); - Array::iterator > normal_traction_it = - normal_traction.begin_reinterpret(nb_quad_facet, spatial_dimension, nb_facet); + // compute average stress on the current quadrature point + Matrix stress_1(facet_stress_it.storage(), + spatial_dimension, + spatial_dimension); - Array::const_iterator sigma_lim_it = sigma_lim.begin(); + Matrix stress_2(facet_stress_it.storage() + sp2, + spatial_dimension, + spatial_dimension); - for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it, ++stress_check_it, - ++normal_traction_it) { - UInt facet = f_filter(f); - if (!facets_check(facet)) continue; + stress_tmp.copy(stress_1); + stress_tmp += stress_2; + stress_tmp /= 2.; - Real mean_stress = std::accumulate(stress_check_it->storage(), - stress_check_it->storage() + nb_quad_facet, - 0.); - mean_stress /= nb_quad_facet; + Vector normal_traction_vec(normal_traction(q)); - if (mean_stress > (*sigma_lim_it - tolerance)) { + // compute normal and effective stress + stress_check(q) = computeEffectiveNorm(stress_tmp, normal, tangent, + normal_traction_vec); + } + + // verify if the effective stress overcomes the threshold + if (stress_check.mean() > (*sigma_lim_it - tolerance)) { f_insertion(facet) = true; + // store the new cohesive material parameters for each quadrature point for (UInt q = 0; q < nb_quad_facet; ++q) { - Real new_sigma = (*stress_check_it)(q); - - Vector ins_s(normal_traction_it->storage() + q * spatial_dimension, - spatial_dimension); + Real new_sigma = stress_check(q); + Vector normal_traction_vec(normal_traction(q)); if (spatial_dimension != 3) - ins_s *= -1.; + normal_traction_vec *= -1.; - sig_c_eff.push_back(new_sigma); - ins_stress.push_back(ins_s); - trac_old.push_back(ins_s); + new_sigmas.push_back(new_sigma); + new_normal_traction.push_back(normal_traction_vec); Real new_delta; // set delta_c in function of G_c or a given delta_c value if (Math::are_float_equal(delta_c, 0.)) new_delta = 2 * G_c / new_sigma; else new_delta = (*sigma_lim_it) / new_sigma * delta_c; - del_c.push_back(new_delta); + new_delta_c.push_back(new_delta); } + } + } - #if defined (AKANTU_DEBUG_TOOLS) && defined(AKANTU_CORE_CXX11) - debug::element_manager.print(debug::_dm_material_cohesive, - [facet, type_facet, nb_quad_facet, &f_stress](const Element & el)->std::string { - std::stringstream sout; - Element facet_el(type_facet, facet, _not_ghost); - if(facet_el == el) { - Array::const_iterator< Vector > stress = f_stress.begin(f_stress.getNbComponent()); - stress += nb_quad_facet * facet; - for (UInt qs = 0; qs < nb_quad_facet; ++qs, ++stress) { - sout << *stress; - } - } - return sout.str(); - }); - #endif + // update material data for the new elements + UInt old_nb_quad_points = sig_c_eff.getSize(); + UInt new_nb_quad_points = new_sigmas.size(); + sig_c_eff.resize(old_nb_quad_points + new_nb_quad_points); + ins_stress.resize(old_nb_quad_points + new_nb_quad_points); + trac_old.resize(old_nb_quad_points + new_nb_quad_points); + del_c.resize(old_nb_quad_points + new_nb_quad_points); + + std::list::iterator new_sig_it = new_sigmas.begin(); + std::list::iterator new_del_it = new_delta_c.begin(); + std::list< Vector >::iterator new_normal_trac_it = new_normal_traction.begin(); + + for (UInt q = 0; q < new_nb_quad_points; ++q, ++new_sig_it, ++new_del_it, ++new_normal_trac_it) { + sig_c_eff(old_nb_quad_points + q) = *new_sig_it; + del_c(old_nb_quad_points + q) = *new_del_it; + for (UInt dim = 0; dim < spatial_dimension; ++dim) { + ins_stress(old_nb_quad_points + q, dim) = (*new_normal_trac_it)(dim); + trac_old(old_nb_quad_points + q, dim) = (*new_normal_trac_it)(dim); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template -inline void MaterialCohesiveLinear::computeEffectiveNorm(const Matrix & stress, +inline Real MaterialCohesiveLinear::computeEffectiveNorm(const Matrix & stress, const Vector & normal, const Vector & tangent, - Vector & normal_traction, - Real & effective_norm) { - AKANTU_DEBUG_IN(); - + Vector & normal_traction) { normal_traction.mul(stress, normal); Real normal_contrib = normal_traction.dot(normal); /// in 3D tangential components must be summed Real tangent_contrib = 0; if (spatial_dimension == 2) { Real tangent_contrib_tmp = normal_traction.dot(tangent); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } else if (spatial_dimension == 3) { for (UInt s = 0; s < spatial_dimension - 1; ++s) { const Vector tangent_v(tangent.storage() + s * spatial_dimension, spatial_dimension); Real tangent_contrib_tmp = normal_traction.dot(tangent_v); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } } tangent_contrib = std::sqrt(tangent_contrib); - normal_contrib = std::max(0., normal_contrib); - effective_norm = std::sqrt(normal_contrib * normal_contrib - + tangent_contrib * tangent_contrib * beta2_inv); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialCohesiveLinear::computeStressNorms(const Array & facet_stress, - Array & stress_check, - Array & normal_traction, - ElementType type_facet) { - AKANTU_DEBUG_IN(); - - Array & facets_check = model->getElementInserter().getCheckFacets(type_facet); - Array & f_filter = facet_filter(type_facet); - - UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbQuadraturePoints(type_facet); - - const Array & tangents = model->getTangents(type_facet); - const Array & normals - = model->getFEEngine("FacetsFEEngine").getNormalsOnQuadPoints(type_facet); - - Real * stress_check_it = stress_check.storage(); - - Array::const_iterator< Vector > normal_begin = - normals.begin(spatial_dimension); - - Array::const_iterator< Vector > tangent_begin = - tangents.begin(tangents.getNbComponent()); - - Array::const_iterator< Matrix > facet_stress_begin = - facet_stress.begin(spatial_dimension, spatial_dimension * 2); - - Array::iterator > normal_traction_it = - normal_traction.begin(spatial_dimension); - - Matrix stress_tmp(spatial_dimension, spatial_dimension); - UInt nb_facet = f_filter.getSize(); - UInt sp2 = spatial_dimension * spatial_dimension; - UInt * current_facet = f_filter.storage(); - - for (UInt f = 0; f < nb_facet; ++f, ++current_facet) { - - if (!facets_check(*current_facet)) { - stress_check_it += nb_quad_facet; - normal_traction_it += nb_quad_facet; - continue; - } - - UInt current_quad = *current_facet * nb_quad_facet; - - for (UInt q = 0; q < nb_quad_facet; ++q, ++stress_check_it, - ++normal_traction_it, ++current_quad) { - - const Vector & normal = normal_begin[current_quad]; - const Vector & tangent = tangent_begin[current_quad]; - const Matrix & facet_stress_it = facet_stress_begin[current_quad]; - - /// compute average stress - Matrix stress_1(facet_stress_it.storage(), - spatial_dimension, - spatial_dimension); - - Matrix stress_2(facet_stress_it.storage() + sp2, - spatial_dimension, - spatial_dimension); - - stress_tmp.copy(stress_1); - stress_tmp += stress_2; - stress_tmp /= 2.; - - /// compute normal, tangential and effective stress - computeEffectiveNorm(stress_tmp, normal, tangent, - *normal_traction_it, *stress_check_it); - } - } - - AKANTU_DEBUG_OUT(); + return std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib * beta2_inv); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::iterator< Vector > traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Array::iterator< Vector > opening_it = opening(el_type, ghost_type).begin(spatial_dimension); Array::iterator< Vector > contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); Array::iterator< Vector > contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::const_iterator< Vector > normal_it = normal.begin(spatial_dimension); Array::iterator< Vector >traction_end = tractions(el_type, ghost_type).end(spatial_dimension); Array::iteratorsigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); Array::iteratordelta_max_it = delta_max(el_type, ghost_type).begin(); Array::iteratordelta_c_it = delta_c_eff(el_type, ghost_type).begin(); Array::iteratordamage_it = damage(el_type, ghost_type).begin(); Array::iterator > insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); Real * memory_space = new Real[2*spatial_dimension]; Vector normal_opening(memory_space, spatial_dimension); 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, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_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 * tangential_opening_norm * beta2_kappa2; bool penetration = normal_opening_norm < -Math::getTolerance(); if (penetration) { /// use penalty coefficient in case of penetration *contact_traction_it = normal_opening; *contact_traction_it *= penalty; *contact_opening_it = normal_opening; /// don't consider penetration contribution for delta *opening_it = tangential_opening; normal_opening.clear(); } else { delta += normal_opening_norm * normal_opening_norm; contact_traction_it->clear(); contact_opening_it->clear(); } delta = std::sqrt(delta); /// update maximum displacement and damage *delta_max_it = std::max(*delta_max_it, delta); *damage_it = std::min(*delta_max_it / *delta_c_it, 1.); /** * 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$ */ if (Math::are_float_equal(*damage_it, 1.)) traction_it->clear(); else if (Math::are_float_equal(*damage_it, 0.)) { if (penetration) traction_it->clear(); else *traction_it = *insertion_stress_it; } else { *traction_it = tangential_opening; *traction_it *= beta2_kappa; *traction_it += normal_opening; AKANTU_DEBUG_ASSERT(*delta_max_it != 0., "Division by zero, tolerance might be too low"); *traction_it *= *sigma_c_it / *delta_max_it * (1. - *damage_it); } } delete [] memory_space; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialCohesiveLinear); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh index 2dfbb16dc..12cbda949 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh @@ -1,165 +1,158 @@ /** * @file material_cohesive_linear.hh * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Tue Jul 29 2014 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_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) * - penalty : stiffness in compression to prevent penetration */ template class MaterialCohesiveLinear : public MaterialCohesive { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialCohesiveLinear(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material parameters virtual void initMaterial(); /// check stress for cohesive elements' insertion virtual void checkInsertion(); protected: /// constitutive law virtual void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost); - /// compute stress norms on quadrature points for each facet for stress check - virtual void computeStressNorms(const Array & facet_stress, - Array & stress_check, - Array & normal_stress, - ElementType type_facet); - /// compute effective stress norm for insertion check - inline void computeEffectiveNorm(const Matrix & stress, + inline Real computeEffectiveNorm(const Matrix & stress, const Vector & normal, const Vector & tangent, - Vector & normal_stress, - Real & effective_norm); + Vector & normal_stress); /** * Scale insertion traction sigma_c according to the volume of the * two elements surrounding a facet * * see the article: F. Zhou and J. F. Molinari "Dynamic crack * propagation with cohesive elements: a methodology to address mesh * dependency" International Journal for Numerical Methods in * Engineering (2004) */ void scaleInsertionTraction(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get sigma_c_eff AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(InsertionTraction, sigma_c_eff, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// beta parameter Real beta; /// beta square inverse to compute effective norm Real beta2_inv; /// mode I fracture energy Real G_c; /// kappa parameter Real kappa; /// constitutive law scalar to compute delta Real beta2_kappa2; /// constitutive law scalar to compute traction Real beta2_kappa; /// penalty coefficient Real penalty; /// reference volume used to scale sigma_c Real volume_s; /// weibull exponent used to scale sigma_c Real m_s; /// critical effective stress RandomInternalField sigma_c_eff; /// effective critical displacement (each element can have a /// different value) CohesiveInternalField delta_c_eff; /// stress at insertion CohesiveInternalField insertion_stress; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_cohesive_linear_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_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 0a2998b6a..d172be9e3 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc @@ -1,618 +1,617 @@ /** * @file material_cohesive.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Jul 29 2014 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" #include "aka_random_generator.hh" #include "shape_cohesive.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model, id), facet_filter("facet_filter", id, this->getMemoryID()), fem_cohesive(&(model.getFEEngineClass("CohesiveFEEngine"))), reversible_energy("reversible_energy", *this), total_energy("total_energy", *this), opening("opening", *this), opening_old("opening (old)", *this), tractions("tractions", *this), tractions_old("tractions (old)", *this), contact_tractions("contact_tractions", *this), contact_opening("contact_opening", *this), delta_max("delta max", *this), use_previous_delta_max(false), damage("damage", *this), - sigma_c("sigma_c", *this) { + sigma_c("sigma_c", *this), + normal(0, spatial_dimension, "normal") { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable, "Critical stress"); this->registerParam("delta_c", delta_c, 0., _pat_parsable | _pat_readable, "Critical displacement"); this->model->getMesh().initElementTypeMapArray(this->element_filter, 1, spatial_dimension, false, _ek_cohesive); if (this->model->getIsExtrinsic()) this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1, spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); this->reversible_energy.initialize(1 ); this->total_energy .initialize(1 ); this->tractions_old .initialize(spatial_dimension); this->tractions .initialize(spatial_dimension); this->opening_old .initialize(spatial_dimension); this->contact_tractions.initialize(spatial_dimension); this->contact_opening .initialize(spatial_dimension); this->opening .initialize(spatial_dimension); this->delta_max .initialize(1 ); this->damage .initialize(1 ); if (model->getIsExtrinsic()) this->sigma_c.initialize(1); if (use_previous_delta_max) delta_max.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Tractions", tractions); #endif Array & 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) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; const Array & shapes = fem_cohesive->getShapes(*it, ghost_type); Array & traction = tractions(*it, ghost_type); Array & contact_traction = contact_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); /// compute @f$t_i N_a@f$ Array * traction_cpy = new Array(nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); Array::iterator > traction_it = traction.begin(spatial_dimension, 1); Array::iterator > contact_traction_it = contact_traction.begin(spatial_dimension, 1); Array::const_iterator > shapes_filtered_begin = shapes.begin(1, size_of_shapes); Array::iterator > traction_cpy_it = traction_cpy->begin(spatial_dimension, size_of_shapes); Matrix traction_tmp(spatial_dimension, 1); for (UInt el = 0; el < nb_element; ++el) { UInt current_quad = elem_filter(el) * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++contact_traction_it, ++current_quad, ++traction_cpy_it) { const Matrix & shapes_filtered = shapes_filtered_begin[current_quad]; traction_tmp.copy(*traction_it); traction_tmp += *contact_traction_it; traction_cpy_it->mul(traction_tmp, 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$ */ Array * int_t_N = new Array(nb_element, spatial_dimension*size_of_shapes, "int_t_N"); fem_cohesive->integrate(*traction_cpy, *int_t_N, spatial_dimension * size_of_shapes, *it, ghost_type, elem_filter); delete traction_cpy; int_t_N->extendComponentsInterlaced(2, int_t_N->getNbComponent()); 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->getFEEngineBoundary().assembleArray(*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(GhostType ghost_type) { AKANTU_DEBUG_IN(); 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 Array & shapes = fem_cohesive->getShapes(*it, ghost_type); Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if(!nb_element) continue; UInt size_of_shapes = shapes.getNbComponent(); Array * shapes_filtered = new Array(nb_element*nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_val = shapes.storage(); Real * shapes_filtered_val = shapes_filtered->storage(); UInt * elem_filter_val = elem_filter.storage(); 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; } /** * 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)); Matrix A(spatial_dimension*size_of_shapes, spatial_dimension*nb_nodes_per_element); for ( UInt i = 0; i < spatial_dimension*size_of_shapes; ++i) { A(i, i) = 1; A(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$ Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix"); // Array * normal = new Array(nb_element * nb_quadrature_points, spatial_dimension, "normal"); - Array normal(nb_quadrature_points, spatial_dimension, "normal"); - - + normal.resize(nb_quadrature_points); computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); tangent_stiffness_matrix->clear(); computeTangentTraction(*it, *tangent_stiffness_matrix, normal, ghost_type); // delete normal; UInt size_at_nt_d_n_a = spatial_dimension*nb_nodes_per_element*spatial_dimension*nb_nodes_per_element; Array * at_nt_d_n_a = new Array (nb_element*nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Array::iterator > shapes_filt_it = shapes_filtered->begin(size_of_shapes); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Array::matrix_iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Array::matrix_iterator At_Nt_D_N_A_end = at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Matrix N (spatial_dimension, spatial_dimension * size_of_shapes); Matrix N_A (spatial_dimension, spatial_dimension * nb_nodes_per_element); Matrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for(; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { N.clear(); /** * store the shapes in voigt notations 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 i = 0; i < spatial_dimension ; ++i) for (UInt n = 0; n < size_of_shapes; ++n) N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n); /** * 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$ **/ N_A.mul(N, A); D_N_A.mul(*D_it, N_A); (*At_Nt_D_N_A_it).mul(D_N_A, N_A); } delete tangent_stiffness_matrix; delete shapes_filtered; Array * K_e = new Array(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive->integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, *it, ghost_type, elem_filter); delete at_nt_d_n_a; model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, *it, ghost_type, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * 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(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Openings", opening); #endif 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) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; UInt nb_quadrature_points = nb_element*fem_cohesive->getNbQuadraturePoints(*it, ghost_type); - Array normal(nb_quadrature_points, spatial_dimension, "normal"); + normal.resize(nb_quadrature_points); /// 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 Array & position, Array & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (type == _cohesive_1d_2) fem_cohesive->computeNormalsOnControlPoints(position, normal, type, ghost_type); else { #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 Array & displacement, Array & 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]; Vector b(memory_space, spatial_dimension); Vector h(memory_space + spatial_dimension, spatial_dimension); for(; it != last_type; ++it) { Array::iterator erev = reversible_energy(*it, _not_ghost).begin(); Array::iterator etot = total_energy(*it, _not_ghost).begin(); Array::vector_iterator traction_it = tractions(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator traction_old_it = tractions_old(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator opening_it = opening(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator opening_old_it = opening_old(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator traction_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; /// update old values it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); GhostType ghost_type = _not_ghost; 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(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate reversible 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 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) { Array 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::getContactEnergy() { AKANTU_DEBUG_IN(); Real econ = 0.; /// integrate contact 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) { Array & el_filter = element_filter(*it, _not_ghost); UInt nb_quad_per_el = fem_cohesive->getNbQuadraturePoints(*it, _not_ghost); UInt nb_quad_points = el_filter.getSize() * nb_quad_per_el; Array contact_energy(nb_quad_points); Array::vector_iterator contact_traction_it = contact_tractions(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator contact_opening_it = contact_opening(*it, _not_ghost).begin(spatial_dimension); /// loop on each quadrature point for (UInt el = 0; el < nb_quad_points; ++contact_traction_it, ++contact_opening_it, ++el) { contact_energy(el) = .5 * contact_traction_it->dot(*contact_opening_it); } econ += fem_cohesive->integrate(contact_energy, *it, _not_ghost, el_filter); } AKANTU_DEBUG_OUT(); return econ; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if (type == "reversible") return getReversibleEnergy(); else if (type == "dissipated") return getDissipatedEnergy(); else if (type == "cohesive contact") return getContactEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ __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 2294ebec7..33489b139 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh @@ -1,246 +1,249 @@ /** * @file material_cohesive.hh * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Jul 29 2014 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" #include "fe_engine_template.hh" #include "aka_common.hh" #include "cohesive_internal_field.hh" #include "cohesive_element_inserter.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 FEEngineTemplate MyFEEngineCohesiveType; public: MaterialCohesive(SolidMechanicsModel& model, const ID & id = ""); virtual ~MaterialCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter virtual void initMaterial(); /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// assemble residual void assembleResidual(GhostType ghost_type = _not_ghost); /// compute reversible and total energies by element void computeEnergies(); /// check stress for cohesive elements' insertion virtual void checkInsertion() { AKANTU_DEBUG_TO_IMPLEMENT(); } /// interpolate stress on given positions for each element (empty /// implemantation to avoid the generic call to be done on cohesive elements) virtual void interpolateStress(__attribute__((unused)) const ElementType type, __attribute__((unused)) Array & result) { }; virtual void computeAllStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) { }; // add the facet to be handled by the material UInt addFacet(const Element & element); protected: virtual void computeTangentTraction(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) const Array & normal, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } void computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type); void computeOpening(const Array & displacement, Array & normal, ElementType type, GhostType ghost_type); template void computeNormal(const Array & position, Array & normal, GhostType ghost_type); /// assemble stiffness void assembleStiffnessMatrix(GhostType ghost_type); /// constitutive law virtual void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// parallelism functions inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* 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); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /// get facet filter AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO(FacetFilter, facet_filter, const ElementTypeMapArray &); // AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// compute contact energy Real getContactEnergy(); /// get energy virtual Real getEnergy(std::string type); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index) { return Material::getEnergy(energy_id, type, index); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// list of facets assigned to this material ElementTypeMapArray facet_filter; /// Link to the cohesive fem object in the model MyFEEngineCohesiveType * fem_cohesive; private: /// reversible energy by quadrature point CohesiveInternalField reversible_energy; /// total energy by quadrature point CohesiveInternalField total_energy; protected: /// opening in all elements and quadrature points CohesiveInternalField opening; /// opening in all elements and quadrature points (previous time step) CohesiveInternalField opening_old; /// traction in all elements and quadrature points CohesiveInternalField tractions; /// traction in all elements and quadrature points (previous time step) CohesiveInternalField tractions_old; /// traction due to contact CohesiveInternalField contact_tractions; /// normal openings for contact tractions CohesiveInternalField contact_opening; /// maximum displacement CohesiveInternalField delta_max; /// tell if the previous delta_max state is needed (in iterative schemes) bool use_previous_delta_max; /// damage CohesiveInternalField damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; /// critical stress RandomInternalField sigma_c; /// critical displacement Real delta_c; + + /// array to temporarily store the normals + Array normal; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_inline_impl.cc" __END_AKANTU__ #include "cohesive_internal_field_tmpl.hh" #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 526a01317..451ab1f95 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc @@ -1,862 +1,868 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Tue May 08 2012 * @date last modification: Fri Sep 05 2014 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include "shape_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "dumpable_inline_impl.hh" #include "material_cohesive.hh" #ifdef AKANTU_USE_IOHELPER # include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass, false, false, true); /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), tangents("tangents", id), stress_on_facet("stress_on_facet", id), facet_stress("facet_stress", id), - facet_material("facet_material", id) { + facet_material("facet_material", id), + cohesive_fe_engine(NULL), + facet_fe_engine(NULL) { AKANTU_DEBUG_IN(); inserter = NULL; #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) facet_synchronizer = NULL; facet_stress_synchronizer = NULL; cohesive_distributed_synchronizer = NULL; global_connectivity = NULL; #endif delete material_selector; material_selector = new DefaultMaterialCohesiveSelector(*this); this->registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, spatial_dimension, _not_ghost, _ek_cohesive); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() { AKANTU_DEBUG_IN(); delete inserter; #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) delete cohesive_distributed_synchronizer; delete facet_synchronizer; delete facet_stress_synchronizer; #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::setTimeStep(Real time_step) { SolidMechanicsModel::setTimeStep(time_step); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFull(const ModelOptions & options) { AKANTU_DEBUG_IN(); const SolidMechanicsModelCohesiveOptions & smmc_options = dynamic_cast(options); this->stress_interpolation = smmc_options.stress_interpolation; this->is_extrinsic = smmc_options.extrinsic; if (!inserter) inserter = new CohesiveElementInserter(mesh, is_extrinsic, synch_parallel, id+":cohesive_element_inserter"); SolidMechanicsModel::initFull(options); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer); #endif if (is_extrinsic) initAutomaticInsertion(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // make sure the material are instantiated if(!are_materials_instantiated) instantiateMaterials(); /// find the first cohesive material UInt 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"); material_selector->setFallback(cohesive_index); // set the facet information in the material in case of dynamic insertion if (is_extrinsic) { const Mesh & mesh_facets = inserter->getMeshFacets(); mesh_facets.initElementTypeMapArray(facet_material, 1, spatial_dimension - 1); Element element; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { element.ghost_type = *gt; Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1, *gt); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, *gt); for(;first != last; ++first) { element.type = *first; Array & f_material = facet_material(*first, *gt); UInt nb_element = mesh_facets.getNbElement(*first, *gt); f_material.resize(nb_element); f_material.set(cohesive_index); for (UInt el = 0; el < nb_element; ++el) { element.element = el; UInt mat_index = (*material_selector)(element); f_material(el) = mat_index; MaterialCohesive & mat = dynamic_cast(*materials[mat_index]); mat.addFacet(element); } } } } else { for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt, _ek_cohesive); Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_cohesive); for(;first != last; ++first) { Array & mat_indexes = this->material_index(*first, *gt); Array & mat_loc_num = this->material_local_numbering(*first, *gt); mat_indexes.set(cohesive_index); mat_loc_num.clear(); } } } SolidMechanicsModel::initMaterials(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModelCohesive::initModel() { AKANTU_DEBUG_IN(); SolidMechanicsModel::initModel(); registerFEEngineObject("CohesiveFEEngine", mesh, spatial_dimension); + cohesive_fe_engine = &getFEEngine("CohesiveFEEngine"); /// add cohesive type connectivity ElementType type = _not_defined; for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType type_ghost = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension, type_ghost); Mesh::type_iterator last = mesh.lastType(spatial_dimension, type_ghost); for (; it != last; ++it) { const Array & connectivity = mesh.getConnectivity(*it, type_ghost); if (connectivity.getSize() != 0) { type = *it; ElementType type_facet = Mesh::getFacetType(type); ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); - getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); - getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); + cohesive_fe_engine->initShapeFunctions(_not_ghost); + cohesive_fe_engine->initShapeFunctions(_ghost); registerFEEngineObject("FacetsFEEngine", mesh.getMeshFacets(), spatial_dimension - 1); + facet_fe_engine = &getFEEngine("FacetsFEEngine"); + if (is_extrinsic) { - getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); - getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); + facet_fe_engine->initShapeFunctions(_not_ghost); + facet_fe_engine->initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::limitInsertion(BC::Axis axis, Real first_limit, Real second_limit) { AKANTU_DEBUG_IN(); inserter->setLimit(axis, first_limit, second_limit); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); inserter->insertIntrinsicElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_stress_synchronizer != NULL) { DataAccessor * data_accessor = this; const ElementTypeMapArray & rank_to_element = synch_parallel->getPrankToElement(); facet_stress_synchronizer->updateFacetStressSynchronizer(*inserter, rank_to_element, *data_accessor); } #endif inserter->getMeshFacets().initElementTypeMapArray(facet_stress, 2 * spatial_dimension * spatial_dimension, spatial_dimension - 1); resizeFacetStress(); /// compute normals on facets computeNormals(); /// allocate stress_on_facet to store element stress on facets mesh.initElementTypeMapArray(stress_on_facet, spatial_dimension * spatial_dimension, spatial_dimension); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_facet_per_elem = Mesh::getNbFacetsPerElement(type); ElementType type_facet = Mesh::getFacetType(type); - UInt nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbQuadraturePoints(type_facet); + UInt nb_quad_per_facet = facet_fe_engine->getNbQuadraturePoints(type_facet); stress_on_facet(type).resize(nb_quad_per_facet * nb_facet_per_elem * nb_element); } if (stress_interpolation) initStressInterpolation(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateAutomaticInsertion() { AKANTU_DEBUG_IN(); inserter->limitCheckFacets(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_stress_synchronizer != NULL) { DataAccessor * data_accessor = this; const ElementTypeMapArray & rank_to_element = synch_parallel->getPrankToElement(); facet_stress_synchronizer->updateFacetStressSynchronizer(*inserter, rank_to_element, *data_accessor); } #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array & position = mesh.getNodes(); ElementTypeMapArray quad_facets("quad_facets", id); mesh_facets.initElementTypeMapArray(quad_facets, spatial_dimension, spatial_dimension - 1); - getFEEngine("FacetsFEEngine").interpolateOnQuadraturePoints(position, quad_facets); + facet_fe_engine->interpolateOnQuadraturePoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray elements_quad_facets("elements_quad_facets", id); mesh.initElementTypeMapArray(elements_quad_facets, spatial_dimension, spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType elem_gt = *gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension, elem_gt); Mesh::type_iterator last = mesh.lastType(spatial_dimension, elem_gt); for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, elem_gt); if (nb_element == 0) continue; /// compute elements' quadrature points and list of facet /// quadrature points positions by element Array & facet_to_element = mesh_facets.getSubelementToElement(type, elem_gt); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); Array & el_q_facet = elements_quad_facets(type, elem_gt); ElementType facet_type = Mesh::getFacetType(type); UInt nb_quad_per_facet = - getFEEngine("FacetsFEEngine").getNbQuadraturePoints(facet_type); + facet_fe_engine->getNbQuadraturePoints(facet_type); el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element global_facet_elem = facet_to_element(el, f); UInt global_facet = global_facet_elem.element; GhostType facet_gt = global_facet_elem.ghost_type; const Array & quad_f = quad_facets(facet_type, facet_gt); 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_f(global_facet * nb_quad_per_facet + q, s); } } } } } } /// loop over non cohesive materials for (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat __attribute__((unused)) = dynamic_cast(*materials[m]); } catch(std::bad_cast&) { /// initialize the interpolation function materials[m]->initElementalFieldInterpolation(elements_quad_facets); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); if (need_initialize) initializeUpdateResidualData(); // f -= fint std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { try { MaterialCohesive & mat = dynamic_cast(**mat_it); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::updateResidual(false); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { try { MaterialCohesive & mat = dynamic_cast(**mat_it); mat.computeEnergies(); } catch (std::bad_cast & bce) { } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = inserter->getMeshFacets(); - getFEEngine("FacetsFEEngine").computeNormalsOnControlPoints(_not_ghost); + facet_fe_engine->computeNormalsOnControlPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = spatial_dimension * (spatial_dimension - 1); mesh_facets.initElementTypeMapArray(tangents, tangent_components, spatial_dimension - 1); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); for (; it != last; ++it) { ElementType facet_type = *it; const Array & normals = - getFEEngine("FacetsFEEngine").getNormalsOnQuadPoints(facet_type); + facet_fe_engine->getNormalsOnQuadPoints(facet_type); UInt nb_quad = normals.getSize(); Array & tang = tangents(facet_type); tang.resize(nb_quad); Real * normal_it = normals.storage(); Real * tangent_it = tang.storage(); /// compute first tangent for (UInt q = 0; q < nb_quad; ++q) { /// if normal is orthogonal to xy plane, arbitrarly define tangent if ( Math::are_float_equal(Math::norm2(normal_it), 0) ) tangent_it[0] = 1; else Math::normal2(normal_it, tangent_it); normal_it += spatial_dimension; tangent_it += tangent_components; } /// compute second tangent (3D case) if (spatial_dimension == 3) { normal_it = normals.storage(); tangent_it = tang.storage(); for (UInt q = 0; q < nb_quad; ++q) { Math::normal3(normal_it, tangent_it, tangent_it + spatial_dimension); normal_it += spatial_dimension; tangent_it += tangent_components; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::averageStressOnFacets(UInt material_index) { AKANTU_DEBUG_IN(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); /// loop over element type for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_quad_per_element = getFEEngine().getNbQuadraturePoints(type); const Array & stress = materials[material_index]->getStress(type); Array & s_on_facet = stress_on_facet(type); UInt nb_facet_per_elem = Mesh::getNbFacetsPerElement(type); ElementType type_facet = Mesh::getFacetType(type); - UInt nb_quad_per_facet = getFEEngine("FacetsFEEngine").getNbQuadraturePoints(type_facet); + UInt nb_quad_per_facet = facet_fe_engine->getNbQuadraturePoints(type_facet); UInt nb_facet_quad_per_elem = nb_quad_per_facet * nb_facet_per_elem; Array::const_iterator > stress_it = stress.begin(spatial_dimension, spatial_dimension); Array::iterator > s_on_facet_it = s_on_facet.begin(spatial_dimension, spatial_dimension); Matrix average_stress(spatial_dimension, spatial_dimension); for (UInt el = 0; el < nb_element; ++el) { /// compute average stress average_stress.clear(); for (UInt q = 0; q < nb_quad_per_element; ++q, ++stress_it) average_stress += *stress_it; average_stress /= nb_quad_per_element; /// store average stress for (UInt q = 0; q < nb_facet_quad_per_elem; ++q, ++s_on_facet_it) *s_on_facet_it = average_stress; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::fillStressOnFacet() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = inserter->getMeshFacets(); UInt sp2 = spatial_dimension * spatial_dimension; UInt sp4 = sp2 * 2; /// loop over materials for (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat __attribute__((unused)) = dynamic_cast(*materials[m]); } catch(std::bad_cast&) { - if (stress_interpolation) + if (stress_interpolation) { /// interpolate stress on facet quadrature points positions - materials[m]->interpolateStress(stress_on_facet); - else + materials[m]->interpolateStressOnFacets(facet_stress); + return; + } else averageStressOnFacets(m); GhostType ghost_type = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last = mesh.lastType(spatial_dimension, ghost_type); /// loop over element type for (; it != last; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) continue; Array & stress_on_f = stress_on_facet(type, ghost_type); /// store the interpolated stresses on the facet_stress vector Array & facet_to_element = mesh_facets.getSubelementToElement(type, ghost_type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); Array::iterator > facet_to_el_it = facet_to_element.begin(nb_facet_per_elem); Array::iterator< Matrix > stress_on_f_it = stress_on_f.begin(spatial_dimension, spatial_dimension); ElementType facet_type = _not_defined; GhostType facet_gt = _casper; Array > * element_to_facet = NULL; Array * f_stress = NULL; Array * f_check = NULL; UInt nb_quad_per_facet = 0; UInt element_rank = 0; Element current_el(type, 0, ghost_type); for (UInt el = 0; el < nb_element; ++el, ++facet_to_el_it) { current_el.element = el; for (UInt f = 0; f < nb_facet_per_elem; ++f) { Element global_facet_elem = (*facet_to_el_it)(f); UInt global_facet = global_facet_elem.element; if (facet_type != global_facet_elem.type || facet_gt != global_facet_elem.ghost_type) { facet_type = global_facet_elem.type; facet_gt = global_facet_elem.ghost_type; element_to_facet = &( mesh_facets.getElementToSubelement(facet_type, facet_gt) ); f_stress = &( facet_stress(facet_type, facet_gt) ); nb_quad_per_facet = - getFEEngine("FacetsFEEngine").getNbQuadraturePoints(facet_type, facet_gt); + facet_fe_engine->getNbQuadraturePoints(facet_type, facet_gt); f_check = &( inserter->getCheckFacets(facet_type, facet_gt) ); } if (!(*f_check)(global_facet)) { stress_on_f_it += nb_quad_per_facet; continue; } for (UInt q = 0; q < nb_quad_per_facet; ++q, ++stress_on_f_it) { element_rank = (*element_to_facet)(global_facet)[0] != current_el; Matrix facet_stress_local(f_stress->storage() + (global_facet * nb_quad_per_facet + q) * sp4 + element_rank * sp2, spatial_dimension, spatial_dimension); facet_stress_local = *stress_on_f_it; } } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::checkCohesiveStress() { AKANTU_DEBUG_IN(); fillStressOnFacet(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_model_cohesive, "Interpolated stresses before", facet_stress); #endif synch_registry->synchronize(_gst_smmc_facets_stress); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_model_cohesive, "Interpolated stresses", facet_stress); #endif for (UInt m = 0; m < materials.size(); ++m) { try { MaterialCohesive & mat_cohesive = dynamic_cast(*materials[m]); /// check which not ghost cohesive elements are to be created mat_cohesive.checkInsertion(); } catch(std::bad_cast&) { } } /// communicate data among processors synch_registry->synchronize(_gst_smmc_facets); /// insert cohesive elements inserter->insertElements(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) updateCohesiveSynchronizers(); #endif SolidMechanicsModel::onElementsAdded(element_list, event); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (cohesive_distributed_synchronizer != NULL) cohesive_distributed_synchronizer->computeAllBufferSizes(*this); #endif /// update shape functions - getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); - getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); + cohesive_fe_engine->initShapeFunctions(_not_ghost); + cohesive_fe_engine->initShapeFunctions(_ghost); if (is_extrinsic) resizeFacetStress(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array & doubled_nodes, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_new_nodes = doubled_nodes.getSize(); Array nodes_list(nb_new_nodes); for (UInt n = 0; n < nb_new_nodes; ++n) nodes_list(n) = doubled_nodes(n, 1); SolidMechanicsModel::onNodesAdded(nodes_list, event); for (UInt n = 0; n < nb_new_nodes; ++n) { UInt old_node = doubled_nodes(n, 0); UInt new_node = doubled_nodes(n, 1); for (UInt dim = 0; dim < spatial_dimension; ++dim) { (*displacement)(new_node, dim) = (*displacement)(old_node, dim); (*velocity) (new_node, dim) = (*velocity) (old_node, dim); (*acceleration)(new_node, dim) = (*acceleration)(old_node, dim); (*blocked_dofs)(new_node, dim) = (*blocked_dofs)(old_node, dim); if (current_position) (*current_position)(new_node, dim) = (*current_position)(old_node, dim); if (increment_acceleration) (*increment_acceleration)(new_node, dim) = (*increment_acceleration)(old_node, dim); if (increment) (*increment)(new_node, dim) = (*increment)(old_node, dim); if (previous_displacement) (*previous_displacement)(new_node, dim) = (*previous_displacement)(old_node, dim); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onEndSolveStep(const AnalysisMethod & method) { AKANTU_DEBUG_IN(); /****************************************************************************** This is required because the Cauchy stress is the stress measure that is used to check the insertion of cohesive elements ******************************************************************************/ std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } 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; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::resizeFacetStress() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = inserter->getMeshFacets(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, ghost_type); for(; it != end; ++it) { ElementType type = *it; UInt nb_facet = mesh_facets.getNbElement(type, ghost_type); UInt nb_quadrature_points = - getFEEngine("FacetsFEEngine").getNbQuadraturePoints(type, ghost_type); + facet_fe_engine->getNbQuadraturePoints(type, ghost_type); UInt new_size = nb_facet * nb_quadrature_points; facet_stress(type, ghost_type).resize(new_size); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump(){ this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ __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 9d6ec740e..b4a3dbc14 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh @@ -1,295 +1,304 @@ /** * @file solid_mechanics_model_cohesive.hh * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Tue Sep 02 2014 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ #include "solid_mechanics_model.hh" #include "solid_mechanics_model_event_handler.hh" #include "cohesive_element_inserter.hh" #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) # include "facet_synchronizer.hh" # include "facet_stress_synchronizer.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelCohesiveOptions : public SolidMechanicsModelOptions { SolidMechanicsModelCohesiveOptions(AnalysisMethod analysis_method = _explicit_lumped_mass, bool extrinsic = false, bool no_init_materials = false, bool stress_interpolation = true) : SolidMechanicsModelOptions(analysis_method, no_init_materials), extrinsic(extrinsic), stress_interpolation(stress_interpolation) {} bool extrinsic; bool stress_interpolation; }; extern const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options; /* -------------------------------------------------------------------------- */ /* Solid Mechanics Model for Cohesive elements */ /* -------------------------------------------------------------------------- */ class SolidMechanicsModelCohesive : public SolidMechanicsModel, public SolidMechanicsModelEventHandler{ /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewCohesiveNodesEvent : public NewNodesEvent { public: AKANTU_GET_MACRO_NOT_CONST(OldNodesList, old_nodes, Array &); AKANTU_GET_MACRO(OldNodesList, old_nodes, const Array &); protected: Array old_nodes; }; typedef FEEngineTemplate MyFEEngineCohesiveType; SolidMechanicsModelCohesive(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model_cohesive", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModelCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step); /// assemble the residual for the explicit scheme virtual void updateResidual(bool need_initialize = true); /// 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(); /// initialize the cohesive model void initFull(const ModelOptions & options = default_solid_mechanics_model_cohesive_options); /// initialize the model void initModel(); /// initialize cohesive material void initMaterials(); /// init facet filters for cohesive materials void initFacetFilter(); /// limit the cohesive element insertion to a given area void limitInsertion(BC::Axis axis, Real first_limit, Real second_limit); /// update automatic insertion after a change in the element inserter void updateAutomaticInsertion(); /// insert intrinsic cohesive elements void insertIntrinsicElements(); private: /// initialize completely the model for extrinsic elements void initAutomaticInsertion(); /// initialize stress interpolation void initStressInterpolation(); /// compute facets' normals void computeNormals(); /// resize facet stress void resizeFacetStress(); /// init facets_check array void initFacetsCheck(); /// fill stress_on_facet void fillStressOnFacet(); /// compute average stress on elements void averageStressOnFacets(UInt material_index); /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded (const Array & nodes_list, const NewNodesEvent & event); virtual void onElementsAdded (const Array & nodes_list, const NewElementsEvent & event); /* ------------------------------------------------------------------------ */ /* SolidMechanicsModelEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: virtual void onEndSolveStep(const AnalysisMethod & method); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get facet mesh AKANTU_GET_MACRO(MeshFacets, mesh.getMeshFacets(), const Mesh &); /// get stress on facets vector AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(StressOnFacets, facet_stress, Real); /// get facet material AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetMaterial, facet_material, UInt); /// get facet material AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetMaterial, facet_material, UInt); /// get facet material AKANTU_GET_MACRO(FacetMaterial, facet_material, const ElementTypeMapArray &); /// @todo THIS HAS TO BE CHANGED AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Tangents, tangents, Real); /// get element inserter AKANTU_GET_MACRO_NOT_CONST(ElementInserter, *inserter, CohesiveElementInserter &); /// get is_extrinsic boolean AKANTU_GET_MACRO(IsExtrinsic, is_extrinsic, bool); + /// get reference to FacetFEEngine + AKANTU_GET_MACRO(FacetFEEngine, *facet_fe_engine, FEEngine &); + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// @todo store tangents when normals are computed: ElementTypeMapArray tangents; /// list of stresses on facet quadrature points for every element ElementTypeMapArray stress_on_facet; /// stress on facets on the two sides by quadrature point ElementTypeMapArray facet_stress; /// material to use if a cohesive element is created on a facet ElementTypeMapArray facet_material; /// stress interpolation flag bool stress_interpolation; bool is_extrinsic; /// cohesive element inserter CohesiveElementInserter * inserter; + /// pointer to the cohesive fe_engine + FEEngine * cohesive_fe_engine; + + /// pointer to the facet fe_engine + FEEngine * facet_fe_engine; + #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) #include "solid_mechanics_model_cohesive_parallel.hh" #endif }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_PARALLEL_COHESIVE_ELEMENT) # include "solid_mechanics_model_cohesive_inline_impl.cc" #endif /* -------------------------------------------------------------------------- */ class DefaultMaterialCohesiveSelector : public DefaultMaterialSelector { public: DefaultMaterialCohesiveSelector(const SolidMechanicsModelCohesive & model) : DefaultMaterialSelector(model.getMaterialByElement()), facet_material(model.getFacetMaterial()), mesh(model.getMesh()) { } inline virtual UInt operator()(const Element & element) { if(Mesh::getKind(element.type) == _ek_cohesive) { try { const Array & cohesive_el_to_facet = mesh.getMeshFacets().getSubelementToElement(element.type, element.ghost_type); bool third_dimension = (mesh.getSpatialDimension() == 3); const Element & facet = cohesive_el_to_facet(element.element, third_dimension); if(facet_material.exists(facet.type, facet.ghost_type)) { return facet_material(facet.type, facet.ghost_type)(facet.element); } else { return MaterialSelector::operator()(element); } } catch (...) { return MaterialSelector::operator()(element); } } else if (Mesh::getSpatialDimension(element.type) == mesh.getSpatialDimension() - 1) { return facet_material(element.type, element.ghost_type)(element.element); } else { return DefaultMaterialSelector::operator()(element); } } private: const ElementTypeMapArray & facet_material; const Mesh & mesh; }; /// 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/test_cohesive_extrinsic.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/test_cohesive_extrinsic.cc index e8b92d2b1..c082b26a0 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,166 +1,166 @@ /** * @file test_cohesive_extrinsic.cc * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Fri Sep 19 2014 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char *argv[]) { initialize("material.dat", argc, argv); debug::setDebugLevel(dblWarning); const UInt spatial_dimension = 2; const UInt max_steps = 1000; Mesh mesh(spatial_dimension); mesh.read("triangle.msh"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); model.limitInsertion(_y, -0.30, -0.20); model.updateAutomaticInsertion(); Real time_step = model.getStableTimeStep()*0.05; model.setTimeStep(time_step); std::cout << "Time step: " << time_step << std::endl; model.assembleMassLumped(); Array & position = mesh.getNodes(); Array & velocity = model.getVelocity(); Array & boundary = model.getBlockedDOFs(); Array & displacement = model.getDisplacement(); // const Array & 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; } /// 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); } model.updateResidual(); model.setBaseName("extrinsic"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity" ); model.addDumpField("acceleration"); model.addDumpField("residual" ); model.addDumpFieldTensor("stress"); model.addDumpField("strain"); model.addDumpField("damage"); model.addDumpFieldToDumper("cohesive elements", "damage"); model.addDumpFieldVectorToDumper("cohesive elements", "displacement"); model.dump(); model.dump("cohesive elements"); // std::ofstream edis("edis.txt"); // std::ofstream erev("erev.txt"); // Array & residual = model.getResidual(); // const Array & 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 < mesh.getNbNodes(); ++n) { if (position(n, 1) > 0.99 || position(n, 1) < -0.99) displacement(n, 1) += disp_update * position(n, 1); } model.checkCohesiveStress(); model.solveStep(); - if(s % 1 == 0) { + if(s % 100 == 0) { model.dump(); model.dump("cohesive elements"); 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(); // mesh.write("mesh_final.msh"); model.dump(); 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::isnan(Ed)) { std::cout << "The dissipated energy is incorrect" << std::endl; finalize(); return EXIT_FAILURE; } finalize(); std::cout << "OK: test_cohesive_extrinsic was passed!" << std::endl; return EXIT_SUCCESS; }