diff --git a/material_laws/damage_law.h b/material_laws/damage_law.h new file mode 100644 index 0000000..78b646c --- /dev/null +++ b/material_laws/damage_law.h @@ -0,0 +1,182 @@ +/* + * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt + * Author: Raphael Prohl + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program 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. + */ + +#ifndef DAMAGE_LAW_H_ +#define DAMAGE_LAW_H_ + +#include "mat_law_interface.h" +#include "lib_disc/function_spaces/grid_function.h" +#include "lib_algebra/cpu_algebra_types.h" + + + +namespace ug{ +namespace SmallStrainMechanics{ + +/// \addtogroup small_strain_mechanics +/// \{ +/// Material Law: Prandtl-Reuss law modelling elastoplastic material behavior where the elastic +/// part is considered as linear. The flow-condition is of von-Mises-type. +/** + * This class implements a material law for small strain elastoplastic material behavior + * + * It is supposed, that the linearized strain tensor could be decomposed additively: + * + * eps = eps_e + eps_p. + * + * The plastic behavior is described by a flow-condition and a flow-rule for the plastic + * evolution (\frac{\partial eps_p){\partial t} = ...). The flow-condition is of + * von-Mises-type and the flow-rule is associative. To treat the plastic equations + * we use the well-established return-mapping-algorithm. Its classical form is valid for the + * 3d-case and the plane strain-case, but not for the plane stress-case! + * + * References: + * + * + * \tparam TDomain + */ + +template +class DamageLaw + : public IMaterialLaw +{ + private: + /// Base class type + typedef IMaterialLaw base_type; + + /// own type + typedef DamageLaw this_type; + + public: + /// World dimension + static const int dim = base_type::dim; + + /// base element type + typedef typename base_type::TBaseElem TBaseElem; + + public: + /// constructor +// DamageLaw(); + DamageLaw( SmartPtr< GridFunction > spPsi0, + SmartPtr< GridFunction > f); + + /// Destructor + ~DamageLaw(){}; + + public: + //////////////////////////// + // INTERFACE-METHODS + //////////////////////////// + void init(); + + /// computes the cauchy stress tensor sigma at an integration point 'ip' + void stressTensor(MathMatrix& stressTens, const size_t ip, + const MathMatrix& GradU); + + /// computes the elasticity tensor; commonly denoted by C + SmartPtr > + elasticityTensor(const size_t ip, const MathMatrix& GradU); + + virtual bool needs_to_add_jac_m(){return false;} + + + virtual void init_internal_vars(TBaseElem* elem, const size_t numIP); + virtual void internal_vars(TBaseElem* elem); + virtual void update_internal_vars(const size_t ip, const MathMatrix& GradU){}; + virtual void attach_internal_vars(typename TDomain::grid_type& grid); + virtual void clear_attachments(typename TDomain::grid_type& grid); + + number& f_on_curr_elem() {return *m_pF_elem;} + number& psi0_on_curr_elem() {return *m_pPsi0_elem;} + + public: + /// set elasticity tensor for orthotropic materials + void set_elasticity_tensor_orthotropic( + const number C11, const number C12, const number C13, + const number C22, const number C23, + const number C33, + const number C44, + const number C55, + const number C66 ); + + /// set hooke elasticity tensor for isotropic materials, (in 2D: plane-strain-case) + void set_hooke_elasticity_tensor(const number lambda, const number mu); + void set_hooke_elasticity_tensor_E_nu(const number E, const number nu); + + public: + using base_type::m_materialConfiguration; + + public: + void strainTensor(MathMatrix& strainTens, const MathMatrix& GradU); + + private: + /// elasticity tensor + SmartPtr > m_spElastTensorFunct; + + + private: + SmartPtr > m_spPsi0; + SmartPtr > m_spF; + + number* m_pF_elem; + number* m_pPsi0_elem; + +/* + // std-vector of InternalVars + struct ElemData{ + number f; + number psi0; + }; + + ElemData* m_pElemData; + + // attachment type: attachment of ElemData + typedef Attachment AElemData; + AElemData m_aElemData; + + // the attachment accessor + typedef Grid::AttachmentAccessor ElemDataAccessor; + ElemDataAccessor m_aaElemData; +*/ +}; + +}// end of namespace SmallStrainMechanics +}// end of namespace ug + +#include "damage_law_impl.h" + +#endif /* DAMAGE_LAW_H_ */ diff --git a/material_laws/damage_law_impl.h b/material_laws/damage_law_impl.h new file mode 100644 index 0000000..b5e7e74 --- /dev/null +++ b/material_laws/damage_law_impl.h @@ -0,0 +1,361 @@ +/* + * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt + * Author: Raphael Prohl + * + * This file is part of UG4. + * + * UG4 is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License version 3 (as published by the + * Free Software Foundation) with the following additional attribution + * requirements (according to LGPL/GPL v3 §7): + * + * (1) The following notice must be displayed in the Appropriate Legal Notices + * of covered and combined works: "Based on UG4 (www.ug4.org/license)". + * + * (2) The following notice must be displayed at a prominent place in the + * terminal output of covered works: "Based on UG4 (www.ug4.org/license)". + * + * (3) The following bibliography is recommended for citation and must be + * preserved in all covered files: + * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively + * parallel geometric multigrid solver on hierarchically distributed grids. + * Computing and visualization in science 16, 4 (2013), 151-164" + * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel + * flexible software system for simulating pde based models on high performance + * computers. Computing and visualization in science 16, 4 (2013), 165-179" + * + * This program 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. + */ + +#ifndef DAMAGE_LAW_IMPL_H_ +#define DAMAGE_LAW_IMPL_H_ + +#include "damage_law.h" + +#include "common/util/string_table_stream.h" + +namespace ug{ +namespace SmallStrainMechanics{ + +/* +template +DamageLaw::DamageLaw(): + IMaterialLaw() +{ + +} +*/ + +template +DamageLaw:: +DamageLaw( SmartPtr > spPsi0, + SmartPtr > f) + : IMaterialLaw(), + m_spPsi0(spPsi0), m_spF(f), + m_pF_elem(NULL), m_pPsi0_elem(NULL) +{ + +} + +template +void +DamageLaw:: +init() +{ + // check, if ElasticityTensor is set + if (m_spElastTensorFunct.invalid()) + UG_THROW("No elasticity tensor set in DamageLaw::init()!"); + + base_type::m_bInit = true; +} + + + +template +void +DamageLaw:: +attach_internal_vars(typename TDomain::grid_type& grid) +{ + /* + grid.template attach_to(m_aElemData); + m_aaElemData.access(grid, m_aElemData); + */ +} + +template +void +DamageLaw:: +clear_attachments(typename TDomain::grid_type& grid) +{ +/* + if(grid.template has_attachment(m_aElemData)){ + grid.template detach_from(m_aElemData); + m_aaElemData.invalidate(); + } +*/ +} + +template +void +DamageLaw:: +init_internal_vars(TBaseElem* elem, const size_t numIP) +{ + +/* + std::vector ind; + const size_t fct = 0; + + if(m_spF->inner_dof_indices(elem, fct, ind) != 1) + UG_THROW("Wrong number dofs"); + + DoFRef(*m_spF, ind[0]) = 1.0; + DoFRef(*m_spPsi0, ind[0]) = 0.0; +*/ +// m_aaElemData[elem].f = 1.0; +// m_aaElemData[elem].psi0 = 0.0; + + if (!base_type::m_bInit){ + base_type::m_bInit = true; + + m_spF->set(1.0); + m_spPsi0->set(0.0); + } +} + +template +inline +void +DamageLaw:: +internal_vars(TBaseElem* elem) +{ + + std::vector ind; + const size_t fct = 0; + + if(m_spF->inner_dof_indices(elem, fct, ind) != 1) + UG_THROW("Wrong number dofs"); + + m_pF_elem = & DoFRef(*m_spF, ind[0]); + m_pPsi0_elem = &DoFRef(*m_spPsi0, ind[0]); + + +// m_pElemData = &m_aaElemData[elem]; +} + + + +template +void +DamageLaw:: +stressTensor(MathMatrix& stressTens, const size_t ip, + const MathMatrix& GradU) +{ + // get linearized strain tensor (eps) at ip + MathMatrix strainTens; + strainTensor(strainTens, GradU); + + // TODO: replace this with general implementation of TensContractMat + for(size_t i = 0; i < (size_t) dim; ++i) + for(size_t j = 0; j < (size_t) dim; ++j) + { + stressTens[i][j] = 0.0; + + for(size_t k = 0; k < (size_t) dim; ++k) + for(size_t l = 0; l < (size_t) dim; ++l) + stressTens[i][j] += (*m_pF_elem) + * (*m_spElastTensorFunct)[i][j][k][l] + * strainTens[k][l]; + + } + + +} + + +template +SmartPtr > +DamageLaw:: +elasticityTensor(const size_t ip, const MathMatrix& GradU_const) +{ + return m_spElastTensorFunct; +} + +template +inline +void +DamageLaw:: +strainTensor(MathMatrix& strainTens, const MathMatrix& GradU) +{ + for(size_t i = 0; i < (size_t) dim; ++i) + for(size_t j = 0; j < (size_t) dim; ++j) + strainTens[i][j] = 0.5 * (GradU[i][j] + GradU[j][i]); +} + +template +void +DamageLaw:: +set_elasticity_tensor_orthotropic( +const number C11, const number C12, const number C13, + const number C22, const number C23, + const number C33, + const number C44, + const number C55, + const number C66 ) +{ + UG_ASSERT( dim==3, "Orthotrope Tensor only for 3 dimensions" ); + + MathTensor4 elastTensorFunct; + + // setze alle wWerte auf 0 + for (size_t i = 0; i < (size_t) dim; ++i) + for (size_t j = 0; j < (size_t) dim; ++j) + for (size_t k = 0; k < (size_t) dim; ++k) + for (size_t l = 0; l < (size_t) dim; ++l) + elastTensorFunct[i][j][k][l] = 0.0; + + + // Tensor mit Werte fuellen + // i j k l + elastTensorFunct[0][0][0][0] = C11; + + elastTensorFunct[0][0][1][1] = C12; + elastTensorFunct[1][1][0][0] = C12; // = C21 + + elastTensorFunct[0][0][2][2] = C13; + elastTensorFunct[2][2][0][0] = C13; // = C31 + + elastTensorFunct[1][1][1][1] = C22; + + elastTensorFunct[1][1][2][2] = C23; + elastTensorFunct[2][2][1][1] = C23; // = C32 + + elastTensorFunct[2][2][2][2] = C33; + + elastTensorFunct[1][2][1][2] = C44; + elastTensorFunct[1][2][2][1] = C44; + + elastTensorFunct[2][1][1][2] = C44; + elastTensorFunct[2][1][2][1] = C44; + + elastTensorFunct[2][0][2][0] = C55; + elastTensorFunct[0][2][2][0] = C55; + elastTensorFunct[2][0][0][2] = C55; + elastTensorFunct[0][2][0][2] = C55; + + elastTensorFunct[0][1][0][1] = C66; + elastTensorFunct[1][0][0][1] = C66; + elastTensorFunct[0][1][1][0] = C66; + elastTensorFunct[1][0][1][0] = C66; + + // remembering the elasticity tensor + SmartPtr > spElastTens(new MathTensor4(elastTensorFunct)); + m_spElastTensorFunct = spElastTens; + + DenseMatrix > mat; + mat(0,0) = C11; mat(0,1) = C12; mat(0,2) = C13; + mat(1,0) = C12; mat(1,1) = C22; mat(1,2) = C23; + mat(2,0) = C13; mat(2,1) = C23; mat(2,2) = C33; + + Invert(mat); + //UG_LOG("S = " mat << "\n"); + number E1 = 1/mat(0,0), E2 = 1/mat(1,1), E3 = 1/mat(2,2); + number v12 = -mat(1, 0)/E1; + //number v21 = -mat(0, 1)/E2; + number v23 = -mat(2, 1)/E2; + //number v32 = -mat(1, 2)/E3; + number v13 = -mat(2, 0)/E1; + //number v31 = -mat(0, 2)/E3; + number G23 = C44; + number G13 = C55; + number G12 = C66; + + StringTableStream sts; + std::stringstream ss; + ss << "Orthrope Elasticity Tensor: \n"; + sts.clear(); + sts << C11 << C12 << C13 << "\n"; + sts << C12 << C22 << C23 << "\n"; + sts << C13 << C23 << C33 << "\n"; + sts << sts.empty_col(3) << C44 << "\n"; + sts << sts.empty_col(4) << C55 << "\n"; + sts << sts.empty_col(5) << C66 << "\n"; + ss << sts; + + ss << "Hooke constants:\n"; + sts.clear(); + sts << " E1 = " << E1 << "E2 = " << E2 << "E3 = " << E3 << "\n"; + sts << " v12 = " << v12 << "v23 = " << v23 << "v13 = " << v13 << "\n"; + sts << " G12 = " << G12 << "G23 = " << G23 << "G13 = " << G13 << "\n"; + ss << sts; + m_materialConfiguration = ss.str(); + UG_LOG("\nset_elasticity_tensor_orthotropic " << ConfigShift(m_materialConfiguration) << "\n"); +} + + +template +void +DamageLaw:: +set_hooke_elasticity_tensor_E_nu(const number E, const number nu) +{ + number lambda = (E*nu) / ((1+nu)*(1-2*nu)); + number mu = E/(2*(1+nu)); + set_hooke_elasticity_tensor(lambda, mu); +} + +template +void +DamageLaw:: +set_hooke_elasticity_tensor(const number lambda, const number mu) +{ + // sets the 'Hooke'-elasticity tensor for isotropic materials + // in 2D this tensor formulation corresponds to the + // plane strain assumption for Hookes`s law + + MathTensor4 elastTensorFunct; + + // filling the constant elasticity tensor + for (size_t i = 0; i < (size_t) dim; ++i) + for (size_t j = 0; j < (size_t) dim; ++j) + for (size_t k = 0; k < (size_t) dim; ++k) + for (size_t l = 0; l < (size_t) dim; ++l) + { + elastTensorFunct[i][j][k][l] = 0.0; + + if ((i == j) && (k == l)) + elastTensorFunct[i][j][k][l] += lambda; + + if ((i == k) && (j == l)) + elastTensorFunct[i][j][k][l] += mu; + + if ((i == l) && (j == k)) + elastTensorFunct[i][j][k][l] += mu; + + } //end (l) + + // remembering the elasticity tensor + SmartPtr > spElastTens(new MathTensor4(elastTensorFunct)); + m_spElastTensorFunct = spElastTens; + + number E = mu * (3.0 * lambda + 2.0 * mu)/(lambda + mu); + number v = 0.5 * lambda/(lambda + mu); + number kappa = lambda + 2.0/3.0 * mu; + + std::stringstream ss; + ss << "Hooke Elasticity Tensor: \n"; + ss << " Lame`s first constant lambda: " << lambda << "\n"; + ss << " Lame`s second constant mue (sometimes 'G', shear modulus) (Schubmodul): " << mu << "\n"; + ss << " This setting equals: \n"; + ss << " young modulus (Elastizitaetsmodul): " << E << "\n"; + ss << " poisson ratio (Querkontraktionszahl) v: " << v << "\n"; + ss << " bulk modulus (Kompressionsmodul): " << kappa << "\n"; + ss << " Elasticity Tensor = " << elastTensorFunct << "\n"; + m_materialConfiguration = ss.str(); +} + +} +} + +#endif /* DAMAGE_LAW_IMPL_H_ */ diff --git a/small_strain_mech.cpp b/small_strain_mech.cpp index 1a47caf..9fb7541 100644 --- a/small_strain_mech.cpp +++ b/small_strain_mech.cpp @@ -31,6 +31,7 @@ */ #include "small_strain_mech.h" +#include "material_laws/damage_law.h" // for various user data #include "bindings/lua/lua_user_data.h" @@ -154,6 +155,15 @@ set_viscous_forces(SmartPtr, dim> > user1, SmartPtr< m_imViscousForces[1].set_data(user2); } + +template +void SmallStrainMechanicsElemDisc:: +set_compress_factor(number val) +{ + m_imCompressIndex.set_data(make_sp(new ConstUserNumber(val))); +} + + ////////////////////////////////// template @@ -206,8 +216,14 @@ init_state_variables(const size_t order) // clears and then attaches the attachments for the internal variables // to the grid SmartPtr dom = this->domain(); + if(dom.invalid()) + UG_THROW("Domain not provided for SmallStrainMechanics::init_state_variables"); + typename TDomain::grid_type& grid = *(dom->grid()); + if(m_spMatLaw.invalid()) + UG_THROW("Material Law not provided for SmallStrainMechanics::init_state_variables"); + m_spMatLaw->clear_attachments(grid); m_spMatLaw->attach_internal_vars(grid); @@ -277,6 +293,7 @@ prep_elem_loop(const ReferenceObjectID roid, const int si) static const int refDim = TElem::dim; m_imVolForce.template set_local_ips(geo.local_ips(), geo.num_ip(), true); m_imDivergence.template set_local_ips(geo.local_ips(), geo.num_ip(), false); + m_imCompressIndex.template set_local_ips(geo.local_ips(), geo.num_ip(), false); m_imViscousForces[0].template set_local_ips(geo.local_ips(), geo.num_ip(), true); m_imViscousForces[1].template set_local_ips(geo.local_ips(), geo.num_ip(), true); } @@ -305,6 +322,7 @@ prep_elem(const LocalVector& u, GridObject* elem, const ReferenceObjectID roid, m_imViscousForces[0].set_global_ips(geo.global_ips(), geo.num_ip()); m_imViscousForces[1].set_global_ips(geo.global_ips(), geo.num_ip()); m_imDivergence.set_global_ips(geo.global_ips(), geo.num_ip()); + m_imCompressIndex.set_global_ips(geo.global_ips(), geo.num_ip()); // set pointer to internal variables of elem @@ -329,6 +347,19 @@ add_jac_A_elem(LocalMatrix& J, const LocalVector& u, m_spElastTensor = m_spMatLaw->elasticityTensor(); } + number f = 1.0; +// number* psi0 = NULL; + DamageLaw* pDamageLaw = dynamic_cast*>(m_spMatLaw.get()); + if(pDamageLaw != NULL) + { + f = pDamageLaw->f_on_curr_elem(); +// psi0 = &pDamageLaw->psi0_on_curr_elem(); + } +// if(psi0) (*psi0) = 0.0; + +// MathMatrix strainTens; +// number volElem = 0.0; + MathMatrix GradU; for (size_t ip = 0; ip < geo.num_ip(); ++ip) { @@ -337,6 +368,9 @@ add_jac_A_elem(LocalMatrix& J, const LocalVector& u, // material law with non constant elasticity tensors m_spMatLaw->template DisplacementGradient(GradU, ip, geo, u); m_spElastTensor = m_spMatLaw->elasticityTensor(ip, GradU); + +// if(pDamageLaw) +// pDamageLaw->strainTensor(strainTens, GradU); } // A) Compute Du:C:Dv = Du:sigma = sigma:Dv @@ -352,9 +386,15 @@ add_jac_A_elem(LocalMatrix& J, const LocalVector& u, for (size_t K = 0; K < (size_t) dim; ++K) for (size_t L = 0; L < (size_t) dim; ++L) { - integrandC += geo.global_grad(ip, a)[K] + integrandC += f + * geo.global_grad(ip, a)[K] * (*m_spElastTensor)[i][K][j][L] * geo.global_grad(ip, b)[L]; + +// if(psi0){ +// (*psi0) += 0.5 * strainTens(i, K) * (*m_spElastTensor)[i][K][j][L] * strainTens(j, L) * geo.weight(ip); +// volElem += geo.weight(ip); +// } } J(i, a, j, b) += integrandC * geo.weight(ip); @@ -364,6 +404,32 @@ add_jac_A_elem(LocalMatrix& J, const LocalVector& u, } //end (j) } } //end(ip) + +// if(psi0){ +// (*psi0) = (*psi0) / volElem; +// } + + if(m_imCompressIndex.data_given()) { + + for(size_t ip = 0; ip < geo.num_ip(); ++ip) + { // for all integration points + + for(size_t i2 = 0; i2 < num_fct(); ++i2){ + for(size_t sh2 = 0; sh2 < geo.num_sh(); ++sh2) + { + + for(size_t sh = 0; sh < geo.num_sh(); ++sh) + { // for all shape functions + for(size_t i = 0; i < num_fct(); ++i) + { // for all components (i=1,2,3) + J(i,sh,i2,sh2) += geo.weight(ip) * m_imCompressIndex[ip] * + geo.global_grad(ip, sh2)[i2] * geo.global_grad(ip, sh)[i]; + } + } + }} + } + } + } // assemble mass jacobian @@ -407,6 +473,27 @@ add_def_A_elem(LocalVector& d, const LocalVector& u, // request geometry const TFEGeom& geo = GeomProvider::get(m_lfeID, m_quadOrder); + + + ////////////// for brittle ////////////////// + number f = 1.0; + number* psi0 = NULL; + DamageLaw* pDamageLaw = dynamic_cast*>(m_spMatLaw.get()); + if(pDamageLaw != NULL) + { + f = pDamageLaw->f_on_curr_elem(); + psi0 = &pDamageLaw->psi0_on_curr_elem(); + } + if(psi0) (*psi0) = 0.0; + + MathMatrix strainTens; + number volElem = 0.0; + ////////////// for brittle ////////////////// + + + + + // a) Cauchy tensor MathMatrix sigma, GradU; for (size_t ip = 0; ip < geo.num_ip(); ++ip) @@ -415,6 +502,13 @@ add_def_A_elem(LocalVector& d, const LocalVector& u, m_spMatLaw->template DisplacementGradient(GradU, ip, geo, u); m_spMatLaw->stressTensor(sigma, ip, GradU); + ////////////// for brittle ////////////////// + m_spElastTensor = m_spMatLaw->elasticityTensor(ip, GradU); + if(pDamageLaw) + pDamageLaw->strainTensor(strainTens, GradU); + ////////////// for brittle ////////////////// + + for (size_t sh = 0; sh < geo.num_sh(); ++sh) {// loop shape functions @@ -428,12 +522,58 @@ add_def_A_elem(LocalVector& d, const LocalVector& u, innerForcesIP += sigma[i][j] * geo.global_grad(ip, sh)[j]; //if (sh>=dim) {UG_LOG("add_def_A_elem: innerForcesIP="<register_import(m_imDivergence); + this->register_import(m_imCompressIndex); this->register_import(m_imVolForce); this->register_import(m_imViscousForces[0]); this->register_import(m_imViscousForces[1]); diff --git a/small_strain_mech.h b/small_strain_mech.h index 6e0fa7a..57f33a9 100644 --- a/small_strain_mech.h +++ b/small_strain_mech.h @@ -167,6 +167,8 @@ class SmallStrainMechanicsElemDisc #endif /// \} + void set_compress_factor(number val); + /** * This methods sets rhs for "viscous stresses" * v0^T (gradPhi +gradPhi^T) v1 @@ -353,6 +355,9 @@ class SmallStrainMechanicsElemDisc /// Data import for the reaction term DataImport m_imDivergence; + /// Data import for the compressibility term + DataImport m_imCompressIndex; + /// data import for viscous forces DataImport, dim > m_imViscousForces[2]; diff --git a/small_strain_mech_plugin.cpp b/small_strain_mech_plugin.cpp index 05244ea..a548445 100644 --- a/small_strain_mech_plugin.cpp +++ b/small_strain_mech_plugin.cpp @@ -34,12 +34,15 @@ #include "bridge/util_domain_algebra_dependent.h" #include "lib_disc/function_spaces/grid_function.h" +#include "lib_disc/common/geometry_util.h" +#include "lib_grid/refinement/refiner_interface.h" #include "small_strain_mech.h" //#include "adaptive_util.h" #include "contact/contact.h" #include "material_laws/hooke.h" +#include "material_laws/damage_law.h" #include "material_laws/prandtl_reuss.h" #include "material_laws/mat_law_interface.h" #include "output_writer/mech_output_writer.h" @@ -51,6 +54,565 @@ using namespace ug::bridge; namespace ug{ namespace SmallStrainMechanics{ + +template +class DamageFunctionUpdater +{ + public: + static const int dim = TDomain::dim; + typedef typename TDomain::grid_type TGrid; + typedef typename grid_dim_traits::element_type TElem; + typedef typename grid_dim_traits::side_type TSide; + + typedef typename TDomain::position_accessor_type TPositionAccessor; + + protected: + void AveragePositions( MathVector& vCenter, + const std::vector >& vCornerCoords) + { + vCenter = vCornerCoords[0]; + for(size_t j = 1; j < vCornerCoords.size(); ++j) + { + vCenter += vCornerCoords[j]; + } + vCenter *= 1./(number)( vCornerCoords.size()); + } + + + void CollectStencilNeighbors(std::vector& vElem, + std::vector& vIndex, + std::vector< MathVector >& vDistance, + TElem* elem, + TGrid& grid, + TPositionAccessor& aaPos, + SmartPtr > spF, + SmartPtr > spPsi0) + { + //static const int numNeighborsToFind = 2*dim + (dim * (dim-1)) / 2; + const size_t fct = 0; + + vElem.clear(); + vIndex.clear(); + vDistance.clear(); + + // get vertices of element + Vertex* const* vVertex = elem->vertices(); + const size_t numVertex = elem->num_vertices(); + + // corner coordinates + std::vector > vCornerCoords; + vCornerCoords.resize(numVertex); + for(size_t vrt = 0; vrt < numVertex; ++vrt){ + vCornerCoords[vrt] = aaPos[ vVertex[vrt] ]; + } + + // element midpoint + MathVector ElemMidPoint; + AveragePositions(ElemMidPoint, vCornerCoords); + +// UG_LOG("############## Element with midpoint " << ElemMidPoint << " ################ \n") + + // get all sides in order + typename TGrid::template traits::secure_container vSide; + grid.associated_elements_sorted(vSide, elem); + + // begin marking + grid.begin_marking(); + + // mark the elem + grid.mark(elem); + + //////////////////////////////////////////////////////////////////////////// + // loop all sides + //////////////////////////////////////////////////////////////////////////// + for(size_t s = 0; s < vSide.size(); ++s){ + +// UG_LOG("######## Boundary Side nr " << s << ": ") + + // get all connected elements + typename TGrid::template traits::secure_container vElemOfSide; + grid.associated_elements(vElemOfSide, vSide[s]); + + // if no elem found: internal error + if(vElemOfSide.size() == 0) + UG_THROW("Huh, why no element?"); + + //////////////////////////////////////////////////////////////////////////// + // introduce mirror element if required + //////////////////////////////////////////////////////////////////////////// + // if only element itself, it must be a boundary element + if(vElemOfSide.size() == 1){ + vElem.push_back(vElemOfSide[0]); + + // add index + std::vector ind; + if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs"); + vIndex.push_back(ind[0]); + + + // add distance + vDistance.resize(vDistance.size()+1); + CollectCornerCoordinates(vCornerCoords, *vSide[s], aaPos); + + MathVector n; + ElementNormal(vSide[s]->reference_object_id(), n, &vCornerCoords[0]); + + ProjectPointToPlane(vDistance.back(), ElemMidPoint, vCornerCoords[0], n); + VecScaleAdd(vDistance.back(), 2.0, ElemMidPoint, -2.0, vDistance.back()); + +// UG_LOG("is boundary, with vDistance: " << vDistance.back() << "\n"); + } + //////////////////////////////////////////////////////////////////////////// + // find all direct face neighbors + //////////////////////////////////////////////////////////////////////////// + // else: add other neighbor and mark + else{ + for(size_t eos = 0; eos < vElemOfSide.size(); ++eos){ + + // neighbor elem + TElem* neighborElem = vElemOfSide[eos]; + + // skip self + if(neighborElem == elem) continue; + + // add neighbor + grid.mark(neighborElem); + vElem.push_back(neighborElem); + + // add index + std::vector ind; + if(spF->inner_dof_indices(neighborElem, fct, ind) != 1) UG_THROW("Wrong number dofs"); + vIndex.push_back(ind[0]); + + // add distance + vDistance.resize(vDistance.size()+1); + CollectCornerCoordinates(vCornerCoords, *neighborElem, aaPos); + AveragePositions(vDistance.back(), vCornerCoords); + VecScaleAdd(vDistance.back(), 1.0, ElemMidPoint, -1.0, vDistance.back()); + +// UG_LOG("has neighbor, c" << vDistance.back() << "\n"); + + } + } + } + + //////////////////////////////////////////////////////////////////////////// + // add additional elements + //////////////////////////////////////////////////////////////////////////// + std::vector vOtherNeighbors; + + + int closest = -1; + number closestDist = std::numeric_limits::max(); + MathVector distance; + +// UG_LOG("Search for vertices-elems: " << numVertex << "\n"); + for(size_t vrt = 0; vrt < numVertex; ++vrt) + { + typename TGrid::template traits::secure_container vVertexNeighbor; + grid.associated_elements(vVertexNeighbor, vVertex[vrt]); + +// UG_LOG(" ++ At vertex "<< vrt << " we have #elems: " << vVertexNeighbor.size() << "\n"); + for(size_t eov = 0; eov < vVertexNeighbor.size(); ++eov) + { + TElem* neighborElem = vVertexNeighbor[eov]; + +// UG_LOG(" ++++ elem "<< eov << " is already marked?: " << grid.is_marked(neighborElem) << "\n"); + if(grid.is_marked(neighborElem)) continue; + + grid.mark(neighborElem); + vOtherNeighbors.push_back(neighborElem); + + + CollectCornerCoordinates(vCornerCoords, *neighborElem, aaPos); + AveragePositions(distance, vCornerCoords); + VecScaleAppend(distance, -1.0, ElemMidPoint); + + number dist = VecTwoNorm(distance); +// UG_LOG(" ++++ ++ dist: " << dist << ", distance: "<< distance << ", ElemMidPoint: "<< ElemMidPoint << "\n") + if(dist < closestDist){ + closest = vOtherNeighbors.size()-1; + } + + } + } + if(dim == 3) UG_THROW("DamageFunctionUpdater: This is 2d only --- extend to 3d by searching for 3 additional neighbors"); + if(closest < 0) UG_THROW("DamageFunctionUpdater: closest not detected.") + + //UG_LOG("vOtherNeighbors.size(): " << vOtherNeighbors.size() << "\n"); + //UG_LOG("closest: " << closest << "\n"); + TElem* otherNeighbor = vOtherNeighbors[closest]; + //UG_LOG("otherNeighbor: " << otherNeighbor << "\n"); + vElem.push_back(otherNeighbor); + + // add index + std::vector ind; + if(spF->inner_dof_indices(otherNeighbor, fct, ind) != 1) UG_THROW("Wrong number dofs"); + vIndex.push_back(ind[0]); + + // add distance + vDistance.resize(vDistance.size()+1); + CollectCornerCoordinates(vCornerCoords, *otherNeighbor, aaPos); + AveragePositions(vDistance.back(), vCornerCoords); + VecScaleAdd(vDistance.back(), 1.0, ElemMidPoint, -1.0, vDistance.back()); + +// UG_LOG("######## Extra Neighbot with vDistance: " << vDistance.back() << "\n") + + // end marking + grid.end_marking(); + } + + + public: + void init( SmartPtr > spF, + SmartPtr > spPsi0) + { + const size_t fct = 0; + + // get domain + SmartPtr domain = spF->domain(); + SmartPtr grid = domain->grid(); + typename TDomain::position_accessor_type& aaPos = domain->position_accessor(); + + // get dof distribution + SmartPtr dd = spF->dd(); + + // get iterators + typename DoFDistribution::traits::iterator iter, iterEnd; + iter = dd->begin(SurfaceView::ALL); // SurfaceView::MG_ALL + iterEnd = dd->end(SurfaceView::ALL); + + // resize result vectors + const size_t numIndex = spF->num_dofs(); + m_vB.resize(numIndex); + m_vDLambda.resize(numIndex); + m_vIndex.resize(numIndex); + + + // storage (for reuse) + std::vector vElem; + std::vector vIndex; + std::vector< MathVector > vDistance; + + // loop all vertices + for(;iter != iterEnd; ++iter) + { + // get vertex + TElem* elem = *iter; + + //////////////////////////////////////////////////////////////////////////// + // Collect suitable neighbors for stencil creation + //////////////////////////////////////////////////////////////////////////// + + CollectStencilNeighbors(vElem, vIndex, vDistance, elem, *grid, aaPos, spF, spPsi0); + + const int numNeighbors = 2*dim + (dim * (dim-1)) / 2; + if(vIndex.size() != numNeighbors || vDistance.size() != numNeighbors) + UG_THROW("Wrong number of neighbors detected: " << vIndex.size()); + + //////////////////////////////////////////////////////////////////////////// + // Create interpolation matrix and invert + //////////////////////////////////////////////////////////////////////////// + + DenseMatrix > BlockInv; + BlockInv.resize(numNeighbors, numNeighbors); + BlockInv = 0.0; + for (size_t j = 0; j < numNeighbors; ++j) + { + for (int d = 0; d < dim; ++d) + BlockInv(j,d) = vDistance[j][d]; + + int cnt = dim; + for (int d1 = 0; d1 < dim-1; ++d1) + for (int d2 = d1+1; d2 < dim; ++d2) + { + BlockInv(j,cnt++) = vDistance[j][d1] * vDistance[j][d2]; + } + + for (int d = 0; d < dim; ++d) + BlockInv(j,cnt++) = 0.5 * vDistance[j][d] * vDistance[j][d]; + } + + if(!Invert(BlockInv)) + UG_THROW("Cannot invert block"); + + //////////////////////////////////////////////////////////////////////////// + // extract second-order derivative subblock + //////////////////////////////////////////////////////////////////////////// + + + // add index + std::vector ind; + if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs"); + const size_t i = ind[0][0]; + + + m_vB[i].resize(dim, numNeighbors); + + m_vDLambda[i] = 0.0; + for (int d = 0; d < dim; ++d){ + for (size_t j = 0; j < numNeighbors; ++j){ + m_vB[i](d,j) = BlockInv( (numNeighbors-dim) + d, j); + + m_vDLambda[i] -= m_vB[i](d,j); + } + } + + + m_vIndex[i].resize(vIndex.size()); + for(size_t k = 0; k < vIndex.size(); ++k) + m_vIndex[i][k] = vIndex[k][0]; + } + } + + + number DLambda(size_t i) {return m_vDLambda[i];} + number Lambda(size_t i, SmartPtr > spF) { + const number fi = (*spF)[i]; + + number res = 0.0; + for (size_t j = 0; j < m_vIndex[i].size(); ++j){ + const number fij = (*spF)[ m_vIndex[i][j] ] - fi; + +// UG_LOG("f(" << i << "," << j << "): " << fij << "\n"); + for (int d = 0; d < dim; ++d){ + res += fij * m_vB[i](d,j); +// UG_LOG("m_vB[" << i << "](" << d << "," << j << "): " << m_vB[i](d,j) << "\n"); + } + } + + return res; + } + + bool solve( SmartPtr > spF, + SmartPtr > spPsi0, + const number beta, const number r, const number eps, const int maxIter) + { + //////////////////////////////////////////////////////////////////////////// + // check if has to be rebuild + //////////////////////////////////////////////////////////////////////////// + + // get approximation space + ConstSmartPtr > approxSpace = spF->approx_space(); + if(approxSpace != spPsi0->approx_space()) + UG_THROW("DamageFunctionUpdater::solve: expected same ApproximationSpace for f and psi0"); + + // check revision counter if grid / approx space has changed since last call + if(m_ApproxSpaceRevision != approxSpace->revision()) + { + // (re-)initialize setting + init(spF, spPsi0); + + // remember revision counter of approx space + m_ApproxSpaceRevision = approxSpace->revision(); + } + + + + //////////////////////////////////////////////////////////////////////////// + // apply newton method + //////////////////////////////////////////////////////////////////////////// + + const size_t numElem = m_vIndex.size(); + + number normPhi = std::numeric_limits::max(); + int iterCnt = 0; + while(normPhi > eps * numElem && (iterCnt++ <= maxIter) ) + { + normPhi = 0.0; + + for(size_t i = 0; i < m_vIndex.size(); ++i) + { + const number lambda = Lambda(i, spF); + number& f = (*spF)[i]; + const number& psi0 = (*spPsi0)[i]; + + number phi = f * ( psi0 - beta * lambda) - r; + + //if( fabs (lambda - 4.0) > 1e-10 ) + // UG_LOG("lambda: " << lambda << "\n"); + (*spPsi0)[i] = lambda; + continue; +// UG_LOG("DLambda: " << DLambda(i) << "\n"); + + + if(phi < eps){ +// phi = 0; +// normPhi += 0.0 * 0.0; + } else { + +/* + UG_LOG(" ############### \n"); + UG_LOG("f : " << f << "\n"); + UG_LOG("psi0 : " << psi0 << "\n"); + UG_LOG("lambda: " << lambda << "\n"); + UG_LOG("DLambda: " << DLambda(i) << "\n"); + UG_LOG("beta : " << beta << "\n"); + UG_LOG("r : " << r << "\n"); + UG_LOG("phi : " << phi << "\n"); + +*/ + f = f - (phi / (psi0 - beta * (lambda + f * DLambda(i)) )); +// UG_LOG ("DamageFunctionUpdater: SETTING new value for f: " << f << "\n"); + + + normPhi += phi*phi; + } + } + + normPhi = sqrt(normPhi); +// UG_LOG(" ######################### (end sweep) ################################### \n"); +// UG_LOG ("DamageFunctionUpdater: normPhi: " << normPhi << "\n"); + } + + if(iterCnt >= maxIter){ + UG_THROW("DamageFunctionUpdater: no convergence after " << iterCnt << " iterations"); + } + + + UG_LOG ("DamageFunctionUpdater: normPhi: " << normPhi << " after " < > > m_vB; + std::vector< number > m_vDLambda; + std::vector< std::vector > m_vIndex; + + + +}; + + + + +template +void MarkDamage( SmartPtr > spF, + SmartPtr > spPsi0, + IRefiner& refiner, + number refineFrac, number coarseFrac, + number avgRefineFactor, number avgCoarsenFactor, + int maxLevel) +{ + static const int dim = TDomain::dim; + typedef typename grid_dim_traits::element_type TElem; + typedef typename DoFDistribution::traits::const_iterator const_iterator; + const int fct = 0; + + /////////////////////////// + // statistic + /////////////////////////// + + number maxValue = std::numeric_limits::min(); + number minValue = std::numeric_limits::max(); + number sumValue = 0.0; + + // loop all elems + const size_t numElem = spF->size(); + for(size_t i = 0; i < numElem; ++i){ + const number val = (*spF)[i] * (*spPsi0)[i]; + + maxValue = std::max(maxValue, val); + minValue = std::min(minValue, val); + sumValue += val; + } + number avgValue = sumValue / numElem; + + + UG_LOG(" +++ numElem: " << numElem << "\n"); + UG_LOG(" +++ maxValue: " << maxValue << "\n"); + UG_LOG(" +++ minValue: " << minValue << "\n"); + UG_LOG(" +++ avgValue: " << avgValue << "\n"); + + + /////////////////////////// + // selection criteria + /////////////////////////// + + number minValueToRefine = std::min( refineFrac * maxValue, avgRefineFactor * avgValue); + number maxValueToCoarsen = std::max( (1 + coarseFrac) * minValue, avgCoarsenFactor * avgValue); + + UG_LOG(" +++ minValueToRefine: " << minValueToRefine << "\n"); + UG_LOG(" +++ maxValueToCoarsen: " << maxValueToCoarsen << "\n"); + + /////////////////////////// + // marking + /////////////////////////// + + // reset counter + int numMarkedRefine = 0, numMarkedCoarse = 0; + + const_iterator iter = spF->template begin(); + const_iterator iterEnd = spF->template end(); + + // loop elements for marking + for(; iter != iterEnd; ++iter) + { + // get element + TElem* elem = *iter; + + std::vector ind; + if(spF->inner_dof_indices(elem, fct, ind) != 1) UG_THROW("Wrong number dofs"); + const number val = DoFRef(*spF, ind[0]) * DoFRef(*spPsi0, ind[0]); + + // marks for refinement + if( val >= minValueToRefine) + if(spF->dd()->multi_grid()->get_level(elem) < maxLevel) + { + refiner.mark(elem, RM_REFINE); + numMarkedRefine++; + } + + // marks for coarsening + if( val <= maxValueToCoarsen) + { + refiner.mark(elem, RM_COARSEN); + numMarkedCoarse++; + } + + } + + /////////////////////////// + // print infos + /////////////////////////// + +#ifdef UG_PARALLEL + if(pcl::NumProcs() > 1) + { + UG_LOG(" +++ Marked for refinement on Proc "< +void HadamardProd( SmartPtr > spFPsi0, + ConstSmartPtr > spF, + ConstSmartPtr > spPsi0) +{ + if(spF->size() != spPsi0->size() || spF->size() != spFPsi0->size()) + UG_THROW("HadamardProd: Size mismatch"); + + for(size_t i = 0; i < spF->size(); ++i){ + (*spFPsi0)[i] = (*spF)[i] * (*spPsi0)[i]; + } +} + + + /** * \defgroup small_strain_mechanics Small Strain Mechanics * \ingroup plugins @@ -148,6 +710,8 @@ static void Domain(Registry& reg, string grp) .add_method("set_div_factor", OVERLOADED_METHOD_PTR(void, T, set_div_factor, (SmartPtr >)), "", "Pressure") .add_method("set_div_factor", OVERLOADED_METHOD_PTR(void, T, set_div_factor, (number)), "", "Pressure") + .add_method("set_compress_factor", OVERLOADED_METHOD_PTR(void, T, set_compress_factor, (number)), "", "Compress-Factor") + .add_method("set_viscous_forces", OVERLOADED_METHOD_PTR(void, T, set_viscous_forces, (SmartPtr, dim> >, SmartPtr, dim> >) ) ,"", "Force field") #ifdef UG_FOR_LUA .add_method("set_volume_forces", OVERLOADED_METHOD_PTR(void, T, set_volume_forces, (const char*)) , "", "Force field") @@ -189,6 +753,41 @@ static void Domain(Registry& reg, string grp) reg.add_class_to_group(name, "HookeLaw", tag); } + + // Damage Law for Linear Elasticity + { + typedef DamageLaw T; + typedef IMaterialLaw TBase; + string name = string("DamageLaw").append(suffix); + reg.add_class_(name, grp) +// .add_constructor() + .template add_constructor >, SmartPtr >)>() + .add_method("set_elasticity_tensor_orthotropic", &T::set_elasticity_tensor_orthotropic, + "", "C11#C12#C13#C22#C23#C33#C44#C55#C66") + //.add_method("set_elasticity_tensor_orthotropic2d", &T::set_elasticity_tensor_orthotropic2d, "", "C11#C12#C22#C33") + .add_method("set_hooke_elasticity_tensor", &T::set_hooke_elasticity_tensor, + "", "lambda#mu") + .add_method("set_hooke_elasticity_tensor_E_nu", &T::set_hooke_elasticity_tensor_E_nu, + "", "E#nu") + .set_construct_as_smart_pointer(true); + reg.add_class_to_group(name, "DamageLaw", tag); + } + + { + typedef DamageFunctionUpdater T; + string name = string("DamageFunctionUpdater").append(suffix); + reg.add_class_(name, grp) + .add_constructor() + .add_method("init", &T::init, "", "") + .add_method("solve", &T::solve, "", "") + .set_construct_as_smart_pointer(true); + reg.add_class_to_group(name, "DamageFunctionUpdater", tag); + + reg.add_function("MarkDamage", &MarkDamage, grp); + reg.add_function("HadamardProd", &HadamardProd, grp); + + } + // Prandtl Reuss Law for small strain ElastoPlasticity { typedef PrandtlReuss T;