diff --git a/packages/damage_non_local.cmake b/packages/damage_non_local.cmake
index 146d0a3a2..830046284 100644
--- a/packages/damage_non_local.cmake
+++ b/packages/damage_non_local.cmake
@@ -1,65 +1,74 @@
#===============================================================================
# @file 25_damage_non_local.cmake
#
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date creation: Fri Jun 15 2012
# @date last modification: Fri Jun 13 2014
#
# @brief package description for non-local materials
#
# @section LICENSE
#
# Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# Akantu is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Akantu. If not, see <http://www.gnu.org/licenses/>.
#
#===============================================================================
package_declare(damage_non_local
DESCRIPTION "Package for Non-local damage constitutives laws Akantu"
DEPENDS lapack)
package_declare_sources(damage_non_local
model/solid_mechanics/materials/material_damage/material_damage_non_local.hh
model/solid_mechanics/materials/material_damage/material_marigo_non_local.hh
model/solid_mechanics/materials/material_damage/material_marigo_non_local_inline_impl.cc
model/solid_mechanics/materials/material_damage/material_mazars_non_local.cc
model/solid_mechanics/materials/material_damage/material_mazars_non_local.hh
model/solid_mechanics/materials/material_non_local.hh
model/solid_mechanics/materials/material_non_local_includes.hh
model/solid_mechanics/materials/material_non_local_inline_impl.cc
- model/solid_mechanics/materials/weight_function.cc
- model/solid_mechanics/materials/weight_function.hh
- model/solid_mechanics/materials/weight_function_tmpl.hh
+ model/solid_mechanics/materials/weight_functions/base_weight_function.hh
+ model/solid_mechanics/materials/weight_functions/base_weight_function_inline_impl.cc
+ model/solid_mechanics/materials/weight_functions/damaged_weight_function.hh
+ model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.cc
+ model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.hh
+ model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.cc
+ model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function.hh
+ model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function_inline_impl.cc
+ model/solid_mechanics/materials/weight_functions/stress_based_weight_function.hh
+ model/solid_mechanics/materials/weight_functions/stress_based_weight_function_tmpl.hh
+ model/solid_mechanics/materials/weight_functions/stress_based_weight_function_inline_impl.cc
+
synchronizer/grid_synchronizer.cc
synchronizer/grid_synchronizer.hh
)
package_declare_material_infos(damage_non_local
LIST AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST
INCLUDE material_non_local_includes.hh
)
package_declare_documentation_files(damage_non_local
manual-constitutive-laws-non_local.tex
manual-appendix-materials-non-local.tex)
package_declare_documentation(damage_non_local
"This package activates the non local damage feature of AKANTU"
"")
diff --git a/src/model/solid_mechanics/materials/material_non_local.hh b/src/model/solid_mechanics/materials/material_non_local.hh
index 905d6f3ff..652920388 100644
--- a/src/model/solid_mechanics/materials/material_non_local.hh
+++ b/src/model/solid_mechanics/materials/material_non_local.hh
@@ -1,185 +1,185 @@
/**
* @file material_non_local.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Wed Aug 31 2011
* @date last modification: Thu Jun 05 2014
*
* @brief Material class that handle the non locality of a law for example damage.
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "material.hh"
#include "aka_grid_dynamic.hh"
#include "fe_engine.hh"
-#include "weight_function.hh"
+#include "base_weight_function.hh"
namespace akantu {
class GridSynchronizer;
}
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_NON_LOCAL_HH__
#define __AKANTU_MATERIAL_NON_LOCAL_HH__
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt dim, template <UInt> class WeightFunction = BaseWeightFunction>
class MaterialNonLocal : public virtual Material {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialNonLocal(SolidMechanicsModel & model, const ID & id = "");
virtual ~MaterialNonLocal();
typedef std::vector< std::pair<QuadraturePoint, QuadraturePoint> > PairList;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize the material computed parameter
virtual void initMaterial();
virtual void updateResidual(GhostType ghost_type);
virtual void computeAllNonLocalStresses(GhostType ghost_type = _not_ghost);
void savePairs(const std::string & filename) const;
void neighbourhoodStatistics(const std::string & filename) const;
protected:
void updatePairList(const ElementTypeMapArray<Real> & quadrature_points_coordinates);
void computeWeights(const ElementTypeMapArray<Real> & quadrature_points_coordinates);
void createCellList(ElementTypeMapArray<Real> & quadrature_points_coordinates);
void cleanupExtraGhostElement(const ElementTypeMap<UInt> & nb_ghost_protected);
void fillCellList(const ElementTypeMapArray<Real> & quadrature_points_coordinates,
const GhostType & ghost_type);
/// constitutive law
virtual void computeNonLocalStresses(GhostType ghost_type = _not_ghost) = 0;
template<typename T>
void weightedAvergageOnNeighbours(const InternalField<T> & to_accumulate,
InternalField<T> & accumulated,
UInt nb_degree_of_freedom,
GhostType ghost_type2 = _not_ghost) const;
virtual inline UInt getNbDataForElements(const Array<Element> & elements,
SynchronizationTag tag) const;
virtual inline void packElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag) const;
virtual inline void unpackElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag);
virtual inline void onElementsAdded(const Array<Element> & element_list);
virtual inline void onElementsRemoved(const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
void registerNonLocalVariable(InternalField<Real> & local,
InternalField<Real> & non_local,
UInt nb_degree_of_freedom) {
ID id = local.getID();
NonLocalVariable & non_local_variable = non_local_variables[id];
non_local_variable.local = &local;
non_local_variable.non_local = &non_local;
non_local_variable.nb_component = nb_degree_of_freedom;
}
AKANTU_GET_MACRO(PairList, pair_list, const PairList &)
Real getRadius() const { return weight_func->getRadius(); }
AKANTU_GET_MACRO(CellList, *spatial_grid, const SpatialGrid<QuadraturePoint> &)
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// the weight function used
WeightFunction<dim> * weight_func;
private:
/**
* the pairs of quadrature points
* 0: not ghost to not ghost
* 1: not ghost to ghost
*/
PairList pair_list[2];
/// the weights associated to the pairs
Array<Real> * pair_weight[2];
/// the regular grid to construct/update the pair lists
SpatialGrid<QuadraturePoint> * spatial_grid;
/// the types of the existing pairs
typedef std::set< std::pair<ElementType, ElementType> > pair_type;
pair_type existing_pairs[2];
/// count the number of calls of computeStress
UInt compute_stress_calls;
struct NonLocalVariable {
InternalField<Real> * local;
InternalField<Real> * non_local;
UInt nb_component;
};
std::map<ID, NonLocalVariable> non_local_variables;
bool is_creating_grid;
GridSynchronizer * grid_synchronizer;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "material_non_local_inline_impl.cc"
__END_AKANTU__
#endif /* __AKANTU_MATERIAL_NON_LOCAL_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_non_local_includes.hh b/src/model/solid_mechanics/materials/material_non_local_includes.hh
index 6dd528ddf..1b9b8d47d 100644
--- a/src/model/solid_mechanics/materials/material_non_local_includes.hh
+++ b/src/model/solid_mechanics/materials/material_non_local_includes.hh
@@ -1,46 +1,51 @@
/**
* @file material_non_local_includes.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Wed Oct 31 2012
* @date last modification: Fri Jun 13 2014
*
* @brief Non local materials includes
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_CMAKE_LIST_MATERIALS
+# include "base_weight_function.hh"
+# include "damaged_weight_function.hh"
+# include "remove_damaged_weight_function.hh"
+# include "remove_damaged_with_damage_rate_weight_function.hh"
+# include "stress_based_weight_function.hh"
# include "material_marigo_non_local.hh"
# include "material_mazars_non_local.hh"
#endif
#define AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST \
((stress_wf, (StressBasedWeightFunction ))) \
((damage_wf, (DamagedWeightFunction ))) \
((remove_wf, (RemoveDamagedWeightFunction))) \
((base_wf, (BaseWeightFunction )))
#define AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST \
((3, (marigo_non_local , MaterialMarigoNonLocal, \
AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST))) \
((2, (mazars_non_local , MaterialMazarsNonLocal )))
diff --git a/src/model/solid_mechanics/materials/weight_function.cc b/src/model/solid_mechanics/materials/weight_function.cc
deleted file mode 100644
index bfe16c82c..000000000
--- a/src/model/solid_mechanics/materials/weight_function.cc
+++ /dev/null
@@ -1,38 +0,0 @@
-/**
- * @file weight_function.cc
- *
- * @author Nicolas Richart <nicolas.richart@epfl.ch>
- *
- * @date creation: Fri Apr 13 2012
- * @date last modification: Tue Nov 06 2012
- *
- * @brief implementation of the weight function classes
- *
- * @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 <http://www.gnu.org/licenses/>.
- *
- */
-
-/* -------------------------------------------------------------------------- */
-
-#include "weight_function.hh"
-
-__BEGIN_AKANTU__
-
-
-__END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/weight_function.hh b/src/model/solid_mechanics/materials/weight_function.hh
deleted file mode 100644
index c66717f53..000000000
--- a/src/model/solid_mechanics/materials/weight_function.hh
+++ /dev/null
@@ -1,374 +0,0 @@
-/**
- * @file weight_function.hh
- *
- * @author Nicolas Richart <nicolas.richart@epfl.ch>
- * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
- *
- * @date creation: Fri Apr 13 2012
- * @date last modification: Thu Jun 05 2014
- *
- * @brief Weight functions for non local materials
- *
- * @section LICENSE
- *
- * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
- * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
- *
- * Akantu is free software: you can redistribute it and/or modify it under the
- * terms of the GNU Lesser General Public License as published by the Free
- * Software Foundation, either version 3 of the License, or (at your option) any
- * later version.
- *
- * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
- * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
- * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
- * details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-/* -------------------------------------------------------------------------- */
-#include "aka_common.hh"
-#include "aka_types.hh"
-#include "solid_mechanics_model.hh"
-#include "parsable.hh"
-#include <cmath>
-#if defined(AKANTU_DEBUG_TOOLS)
-#include "aka_debug_tools.hh"
-#include <string>
-#endif
-
-#include "material_list.hh"
-
-/* -------------------------------------------------------------------------- */
-#include <vector>
-#include "material_damage.hh"
-
-#ifndef __AKANTU_WEIGHT_FUNCTION_HH__
-#define __AKANTU_WEIGHT_FUNCTION_HH__
-
-__BEGIN_AKANTU__
-
-/* -------------------------------------------------------------------------- */
-/* Normal weight function */
-/* -------------------------------------------------------------------------- */
-template<UInt spatial_dimension>
-class BaseWeightFunction : public Parsable {
-public:
- BaseWeightFunction(Material & material, const std::string & type = "base") :
- Parsable(_st_non_local, "weight_function:" + type), material(material), type(type) {
- this->registerParam("radius" , R , 100.,
- _pat_parsable | _pat_readable , "Non local radius");
- this->registerParam("update_rate" , update_rate, 0U ,
- _pat_parsmod, "Update frequency");
- }
-
- virtual ~BaseWeightFunction() {}
-
- virtual void init() { R2 = R * R; };
-
- virtual void updateInternals(__attribute__((unused)) const ElementTypeMapArray<Real> & quadrature_points_coordinates) {};
-
- /* ------------------------------------------------------------------------ */
- inline void setRadius(Real radius) { R = radius; R2 = R * R; }
-
- /* ------------------------------------------------------------------------ */
- inline void selectType(__attribute__((unused)) ElementType type1,
- __attribute__((unused)) GhostType ghost_type1,
- __attribute__((unused)) ElementType type2,
- __attribute__((unused)) GhostType ghost_type2) {
- }
-
- /* ------------------------------------------------------------------------ */
- inline Real operator()(Real r,
- const __attribute__((unused)) QuadraturePoint & q1,
- const __attribute__((unused)) QuadraturePoint & q2) {
- Real w = 0;
- if(r <= R) {
- Real alpha = (1. - r*r / R2);
- w = alpha * alpha;
- // *weight = 1 - sqrt(r / radius);
- }
- return w;
- }
-
- void printself(std::ostream & stream, int indent) const {
- std::string space;
- for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
- stream << space << "WeightFunction " << type << " [" << std::endl;
- Parsable::printself(stream, indent);
- stream << space << "]" << std::endl;
- }
-
-public:
- Real getRadius() { return R; }
- UInt getUpdateRate() { return update_rate; }
-
-public:
- virtual UInt getNbDataForElements(__attribute__((unused)) const Array<Element> & elements,
- __attribute__((unused)) SynchronizationTag tag) const {
- return 0;
- }
-
- virtual inline void packElementData(__attribute__((unused)) CommunicationBuffer & buffer,
- __attribute__((unused)) const Array<Element> & elements,
- __attribute__((unused)) SynchronizationTag tag) const {}
-
- virtual inline void unpackElementData(__attribute__((unused)) CommunicationBuffer & buffer,
- __attribute__((unused)) const Array<Element> & elements,
- __attribute__((unused)) SynchronizationTag tag) {}
-protected:
- Material & material;
-
- Real R;
- Real R2;
-
- UInt update_rate;
-
- const std::string type;
-};
-
-/* -------------------------------------------------------------------------- */
-/* Damage weight function */
-/* -------------------------------------------------------------------------- */
-template<UInt spatial_dimension>
-class DamagedWeightFunction : public BaseWeightFunction<spatial_dimension> {
-public:
- DamagedWeightFunction(Material & material) : BaseWeightFunction<spatial_dimension>(material, "damaged") {
- //AKANTU_DEBUG_ASSERT(dynamic_cast<MaterialDamage<spatial_dimension> *>(&material) != NULL, "This weight function works only with damage materials!");
- }
-
- inline void selectType(__attribute__((unused)) ElementType type1,
- __attribute__((unused)) GhostType ghost_type1,
- ElementType type2,
- GhostType ghost_type2) {
- selected_damage = &(this->material.template getArray<Real>("damage", type2, ghost_type2));
- }
-
- inline Real operator()(Real r,
- const __attribute__((unused)) QuadraturePoint & q1,
- const QuadraturePoint & q2) {
- UInt quad = q2.global_num;
- Real D = (*selected_damage)(quad);
- Real Radius_t = 0;
- Real Radius_init = this->R2;
-
-// if(D <= 0.5)
-// {
-// Radius_t = 2*D*Radius_init;
-// }
-// else
-// {
-// Radius_t = 2*Radius_init*(1-D);
-// }
-//
-
- Radius_t = Radius_init*(1-D);
-
-
- Radius_init *= Radius_init;
- Radius_t *= Radius_t;
-
- if(Radius_t < Math::getTolerance()) {
- Radius_t = 0.001*Radius_init;
- }
-
- Real expb = (2*std::log(0.51))/(std::log(1.0-0.49*Radius_t/Radius_init));
- Int expb_floor=std::floor(expb);
- Real b = expb_floor + expb_floor%2;
- Real alpha = std::max(0., 1. - r*r / Radius_init);
- Real w = std::pow(alpha,b);
- return w;
- }
-
-private:
- const Array<Real> * selected_damage;
-};
-
-/* -------------------------------------------------------------------------- */
-/* Remove damaged weight function */
-/* -------------------------------------------------------------------------- */
-template<UInt spatial_dimension>
-class RemoveDamagedWeightFunction : public BaseWeightFunction<spatial_dimension> {
-public:
- RemoveDamagedWeightFunction(Material & material) : BaseWeightFunction<spatial_dimension>(material, "remove_damaged") {
- this->registerParam("damage_limit", this->damage_limit, 1., _pat_parsable, "Damage Threshold");
- }
-
- inline void selectType(__attribute__((unused)) ElementType type1,
- __attribute__((unused)) GhostType ghost_type1,
- ElementType type2,
- GhostType ghost_type2) {
- selected_damage = &(this->material.template getArray<Real>("damage", type2, ghost_type2));
- // selected_damage = &mat.getDamage(type2, ghost_type2);
- }
-
- inline Real operator()(Real r,
- const __attribute__((unused)) QuadraturePoint & q1,
- const QuadraturePoint & q2) {
- UInt quad = q2.global_num;
-
- if(q1 == q2) return 1.;
-
- Real D = (*selected_damage)(quad);
- Real w = 0.;
- if(D < damage_limit) {
- Real alpha = std::max(0., 1. - r*r / this->R2);
- w = alpha * alpha;
- }
- return w;
- }
-
-
- virtual UInt getNbDataForElements(const Array<Element> & elements,
- SynchronizationTag tag) const {
- if(tag == _gst_mnl_weight)
- return this->material.getModel().getNbQuadraturePoints(elements) * sizeof(Real);
-
- return 0;
- }
-
- virtual inline void packElementData(CommunicationBuffer & buffer,
- const Array<Element> & elements,
- SynchronizationTag tag) const {
- if(tag == _gst_mnl_weight) {
- ElementTypeMapArray<Real> & damage = this->material.template getInternal<Real>("damage");
- this->material.packElementDataHelper(damage,
- buffer,
- elements);
- }
- }
-
- virtual inline void unpackElementData(CommunicationBuffer & buffer,
- const Array<Element> & elements,
- SynchronizationTag tag) {
- if(tag == _gst_mnl_weight) {
- ElementTypeMapArray<Real> & damage = this->material.template getInternal<Real>("damage");
- this->material.unpackElementDataHelper(damage,
- buffer,
- elements);
- }
- }
-
-
-private:
- /// limit at which a point is considered as complitely broken
- Real damage_limit;
-
- /// internal pointer to the current damage vector
- const Array<Real> * selected_damage;
-};
-/* -------------------------------------------------------------------------- */
-/* Remove damaged with damage rate weight function */
-/* -------------------------------------------------------------------------- */
-
-template<UInt spatial_dimension>
-class RemoveDamagedWithDamageRateWeightFunction : public BaseWeightFunction<spatial_dimension> {
-public:
- RemoveDamagedWithDamageRateWeightFunction(Material & material) : BaseWeightFunction<spatial_dimension>(material, "remove_damage_with_damage_rate") {
- this->registerParam("damage_limit", this->damage_limit_with_damage_rate, 1, _pat_parsable, "Damage Threshold");
- }
-
- inline void selectType(__attribute__((unused)) ElementType type1,
- __attribute__((unused)) GhostType ghost_type1,
- ElementType type2,
- GhostType ghost_type2) {
- selected_damage_with_damage_rate = &(this->material.template getArray<Real>("damage",type2, ghost_type2));
- selected_damage_rate_with_damage_rate = &(this->material.template getArray<Real>("damage-rate",type2, ghost_type2));
- }
-
- inline Real operator()(Real r,
- const __attribute__((unused)) QuadraturePoint & q1,
- const QuadraturePoint & q2) {
- UInt quad = q2.global_num;
-
- if(q1.global_num == quad) return 1.;
-
- Real D = (*selected_damage_with_damage_rate)(quad);
- Real w = 0.;
- Real alphaexp = 1.;
- Real betaexp = 2.;
- if(D < damage_limit_with_damage_rate) {
- Real alpha = std::max(0., 1. - pow((r*r / this->R2),alphaexp));
- w = pow(alpha, betaexp);
- }
-
- return w;
- }
-
-private:
- /// limit at which a point is considered as complitely broken
- Real damage_limit_with_damage_rate;
-
- /// internal pointer to the current damage vector
- const Array<Real> * selected_damage_with_damage_rate;
-
- /// internal pointer to the current damage rate vector
- const Array<Real> * selected_damage_rate_with_damage_rate;
-};
-
-/* -------------------------------------------------------------------------- */
-/* Stress Based Weight */
-/* -------------------------------------------------------------------------- */
-template<UInt spatial_dimension>
-class StressBasedWeightFunction : public BaseWeightFunction<spatial_dimension> {
-public:
- StressBasedWeightFunction(Material & material);
-
- void init();
-
- virtual void updateInternals(__attribute__((unused)) const ElementTypeMapArray<Real> & quadrature_points_coordinates) {
- updatePrincipalStress(_not_ghost);
- updatePrincipalStress(_ghost);
- };
-
- void updatePrincipalStress(GhostType ghost_type);
-
- inline void updateQuadraturePointsCoordinates(ElementTypeMapArray<Real> & quadrature_points_coordinates);
-
- inline void selectType(ElementType type1, GhostType ghost_type1,
- ElementType type2, GhostType ghost_type2);
-
- inline Real operator()(Real r,
- const QuadraturePoint & q1,
- const QuadraturePoint & q2);
-
- inline Real computeRhoSquare(Real r,
- Vector<Real> & eigs,
- Matrix<Real> & eigenvects,
- Vector<Real> & x_s);
-
-private:
- Real ft;
-
- InternalField<Real> stress_diag;
- Array<Real> * selected_stress_diag;
- InternalField<Real> stress_base;
- Array<Real> * selected_stress_base;
-
-
- // InternalField<Real> quadrature_points_coordinates;
- Array<Real> * selected_position_1;
- Array<Real> * selected_position_2;
-
- InternalField<Real> characteristic_size;
- Array<Real> * selected_characteristic_size;
-};
-
-
-template<UInt spatial_dimension>
-inline std::ostream & operator <<(std::ostream & stream,
- const BaseWeightFunction<spatial_dimension> & _this)
-{
- _this.printself(stream);
- return stream;
-}
-
-
-#include "weight_function_tmpl.hh"
-
-__END_AKANTU__
-
-#endif /* __AKANTU_WEIGHT_FUNCTION_HH__ */
diff --git a/src/model/solid_mechanics/materials/weight_functions/base_weight_function.hh b/src/model/solid_mechanics/materials/weight_functions/base_weight_function.hh
new file mode 100644
index 000000000..046b66516
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/base_weight_function.hh
@@ -0,0 +1,177 @@
+/**
+ * @file base_weight_function.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief Base weight function for non local materials
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * Akantu is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "aka_common.hh"
+#include "aka_types.hh"
+#include "solid_mechanics_model.hh"
+#include "parsable.hh"
+#include <cmath>
+#if defined(AKANTU_DEBUG_TOOLS)
+#include "aka_debug_tools.hh"
+#include <string>
+#endif
+
+#include "material_list.hh"
+
+/* -------------------------------------------------------------------------- */
+#include <vector>
+#include "material_damage.hh"
+
+#ifndef __AKANTU_BASE_WEIGHT_FUNCTION_HH__
+#define __AKANTU_BASE_WEIGHT_FUNCTION_HH__
+
+__BEGIN_AKANTU__
+
+/* -------------------------------------------------------------------------- */
+/* Normal weight function */
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+class BaseWeightFunction : public Parsable {
+public:
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+ BaseWeightFunction(Material & material, const std::string & type = "base") :
+ Parsable(_st_non_local, "weight_function:" + type), material(material), type(type) {
+ this->registerParam("radius" , R , 100.,
+ _pat_parsable | _pat_readable , "Non local radius");
+ this->registerParam("update_rate" , update_rate, 0U ,
+ _pat_parsmod, "Update frequency");
+ }
+
+ virtual ~BaseWeightFunction() {}
+
+ /* -------------------------------------------------------------------------- */
+ /* Methods */
+ /* -------------------------------------------------------------------------- */
+ /// initialize the weight function
+ virtual inline void init();
+
+ /// update the internal parameters
+ virtual void updateInternals(__attribute__((unused)) const ElementTypeMapArray<Real> & quadrature_points_coordinates) {};
+
+ /* ------------------------------------------------------------------------ */
+ /// set the non-local radius
+ inline void setRadius(Real radius);
+
+ /* ------------------------------------------------------------------------ */
+ /// function for optimization to preselect types for the
+ /// ElementTypeMapArrays (this function is currently not used)
+ inline void selectType(__attribute__((unused)) ElementType type1,
+ __attribute__((unused)) GhostType ghost_type1,
+ __attribute__((unused)) ElementType type2,
+ __attribute__((unused)) GhostType ghost_type2) {
+ }
+
+ /* ------------------------------------------------------------------------ */
+ /// compute the weight for a given distance between two quadrature points
+ inline Real operator()(Real r,
+ const __attribute__((unused)) QuadraturePoint & q1,
+ const __attribute__((unused)) QuadraturePoint & q2);
+
+ /// print function
+ void printself(std::ostream & stream, int indent) const {
+ std::string space;
+ for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
+ stream << space << "WeightFunction " << type << " [" << std::endl;
+ Parsable::printself(stream, indent);
+ stream << space << "]" << std::endl;
+ }
+
+ /* -------------------------------------------------------------------------- */
+ /* Accessors */
+ /* -------------------------------------------------------------------------- */
+
+public:
+ /// get the radius
+ Real getRadius() { return R; }
+ /// get the update rate
+ UInt getUpdateRate() { return update_rate; }
+
+public:
+
+ /* ------------------------------------------------------------------------ */
+ /* Data Accessor inherited members */
+ /* ------------------------------------------------------------------------ */
+
+ virtual UInt getNbDataForElements(__attribute__((unused)) const Array<Element> & elements,
+ __attribute__((unused)) SynchronizationTag tag) const {
+ return 0;
+ }
+
+ virtual inline void packElementData(__attribute__((unused)) CommunicationBuffer & buffer,
+ __attribute__((unused)) const Array<Element> & elements,
+ __attribute__((unused)) SynchronizationTag tag) const {}
+
+ virtual inline void unpackElementData(__attribute__((unused)) CommunicationBuffer & buffer,
+ __attribute__((unused)) const Array<Element> & elements,
+ __attribute__((unused)) SynchronizationTag tag) {}
+protected:
+ /* ------------------------------------------------------------------------ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------ */
+ /// the material for which the weights are computed
+ Material & material;
+
+ /// the non-local radius
+ Real R;
+
+ /// the non-local radius squared
+ Real R2;
+
+ /// the update rate
+ UInt update_rate;
+
+ /// name of the type of weight function
+ const std::string type;
+};
+
+
+template<UInt spatial_dimension>
+inline std::ostream & operator <<(std::ostream & stream,
+ const BaseWeightFunction<spatial_dimension> & _this)
+{
+ _this.printself(stream);
+ return stream;
+}
+
+
+#if defined (AKANTU_INCLUDE_INLINE_IMPL)
+# include "base_weight_function_inline_impl.cc"
+#endif
+
+__END_AKANTU__
+
+/* -------------------------------------------------------------------------- */
+
+
+#endif /* __AKANTU_BASE_WEIGHT_FUNCTION_HH__ */
diff --git a/src/model/solid_mechanics/materials/weight_functions/base_weight_function_inline_impl.cc b/src/model/solid_mechanics/materials/weight_functions/base_weight_function_inline_impl.cc
new file mode 100644
index 000000000..a37a92419
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/base_weight_function_inline_impl.cc
@@ -0,0 +1,61 @@
+/**
+ * @file base_weight_function_inline_impl.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief Implementation of inline function of base weight function
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline void BaseWeightFunction<spatial_dimension>::init() {
+ /// compute R^2 for a given non-local radius
+ this->R2 = this->R *this-> R;
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline void BaseWeightFunction<spatial_dimension>::setRadius(Real radius) {
+ /// set the non-local radius and update R^2 accordingly
+ this->R = radius;
+ this->R2 = this->R * this->R;
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline Real BaseWeightFunction<spatial_dimension>:: operator()(Real r,
+ const __attribute__((unused)) QuadraturePoint & q1,
+ const __attribute__((unused)) QuadraturePoint & q2) {
+ /// initialize the weight
+ Real w = 0;
+ /// compute weight for given r
+ if(r <= this->R) {
+ Real alpha = (1. - r*r / this->R2);
+ w = alpha * alpha;
+ // *weight = 1 - sqrt(r / radius);
+ }
+ return w;
+}
diff --git a/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function.hh b/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function.hh
new file mode 100644
index 000000000..87dbcc153
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function.hh
@@ -0,0 +1,81 @@
+/**
+ * @file damaged_weight_function.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief Damaged weight function for non local materials
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * Akantu is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "base_weight_function.hh"
+/* -------------------------------------------------------------------------- */
+#ifndef __AKANTU_DAMAGED_WEIGHT_FUNCTION_HH__
+#define __AKANTU_DAMAGED_WEIGHT_FUNCTION_HH__
+
+__BEGIN_AKANTU__
+
+/* -------------------------------------------------------------------------- */
+/* Damage weight function */
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+class DamagedWeightFunction : public BaseWeightFunction<spatial_dimension> {
+public:
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+ DamagedWeightFunction(Material & material) : BaseWeightFunction<spatial_dimension>(material, "damaged") {
+ //AKANTU_DEBUG_ASSERT(dynamic_cast<MaterialDamage<spatial_dimension> *>(&material) != NULL, "This weight function works only with damage materials!");
+ }
+
+ /* -------------------------------------------------------------------------- */
+ /* Base Weight Function inherited methods */
+ /* -------------------------------------------------------------------------- */
+ inline void selectType(__attribute__((unused)) ElementType type1,
+ __attribute__((unused)) GhostType ghost_type1,
+ ElementType type2,
+ GhostType ghost_type2);
+
+ inline Real operator()(Real r,
+ const __attribute__((unused)) QuadraturePoint & q1,
+ const QuadraturePoint & q2);
+
+private:
+
+ /* ------------------------------------------------------------------------ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------ */
+ /// member for optimization when types for ElementTypeMaps are preselected (currently not in use)
+ const Array<Real> * selected_damage;
+};
+
+
+#if defined (AKANTU_INCLUDE_INLINE_IMPL)
+# include "damaged_weight_function_inline_impl.cc"
+#endif
+
+__END_AKANTU__
+
+#endif /* __AKANTU_DAMAGED_WEIGHT_FUNCTION_HH__ */
diff --git a/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.cc b/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.cc
new file mode 100644
index 000000000..d535e11b8
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.cc
@@ -0,0 +1,79 @@
+/**
+ * @file damaged_weight_function_inline_impl.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief Implementation of inline function of damaged weight function
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline void DamagedWeightFunction<spatial_dimension>::selectType(__attribute__((unused)) ElementType type1,
+ __attribute__((unused)) GhostType ghost_type1,
+ ElementType type2,
+ GhostType ghost_type2) {
+ /// preselect the damage array for a given type: for optimization
+ selected_damage = &(this->material.template getArray<Real>("damage", type2, ghost_type2));
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline Real DamagedWeightFunction<spatial_dimension>::operator()(Real r,
+ const __attribute__((unused)) QuadraturePoint & q1,
+ const QuadraturePoint & q2) {
+ /// compute the weight
+ UInt quad = q2.global_num;
+ Real D = (*selected_damage)(quad);
+ Real Radius_t = 0;
+ Real Radius_init = this->R2;
+
+ // if(D <= 0.5)
+ // {
+ // Radius_t = 2*D*Radius_init;
+ // }
+ // else
+ // {
+ // Radius_t = 2*Radius_init*(1-D);
+ // }
+ //
+
+ Radius_t = Radius_init*(1-D);
+
+
+ Radius_init *= Radius_init;
+ Radius_t *= Radius_t;
+
+ if(Radius_t < Math::getTolerance()) {
+ Radius_t = 0.001*Radius_init;
+ }
+
+ Real expb = (2*std::log(0.51))/(std::log(1.0-0.49*Radius_t/Radius_init));
+ Int expb_floor=std::floor(expb);
+ Real b = expb_floor + expb_floor%2;
+ Real alpha = std::max(0., 1. - r*r / Radius_init);
+ Real w = std::pow(alpha,b);
+ return w;
+}
diff --git a/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.hh b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.hh
new file mode 100644
index 000000000..aab427eb1
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.hh
@@ -0,0 +1,100 @@
+
+/**
+ * @file remove_damaged_weight_function.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief Removed damaged weight function for non local materials
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * Akantu is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "base_weight_function.hh"
+/* -------------------------------------------------------------------------- */
+#ifndef __AKANTU_REMOVE_DAMAGED_WEIGHT_FUNCTION_HH__
+#define __AKANTU_REMOVE_DAMAGED_WEIGHT_FUNCTION_HH__
+
+__BEGIN_AKANTU__
+
+/* -------------------------------------------------------------------------- */
+/* Remove damaged weight function */
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+class RemoveDamagedWeightFunction : public BaseWeightFunction<spatial_dimension> {
+public:
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+ RemoveDamagedWeightFunction(Material & material) : BaseWeightFunction<spatial_dimension>(material, "remove_damaged") {
+ this->registerParam("damage_limit", this->damage_limit, 1., _pat_parsable, "Damage Threshold");
+ }
+
+ /* -------------------------------------------------------------------------- */
+ /* Base Weight Function inherited methods */
+ /* -------------------------------------------------------------------------- */
+ inline void selectType(__attribute__((unused)) ElementType type1,
+ __attribute__((unused)) GhostType ghost_type1,
+ ElementType type2,
+ GhostType ghost_type2);
+
+ inline Real operator()(Real r,
+ const __attribute__((unused)) QuadraturePoint & q1,
+ const QuadraturePoint & q2);
+
+ /* ------------------------------------------------------------------------ */
+ /* Data Accessor inherited members */
+ /* ------------------------------------------------------------------------ */
+
+ virtual inline UInt getNbDataForElements(const Array<Element> & elements,
+ SynchronizationTag tag) const;
+
+ virtual inline void packElementData(CommunicationBuffer & buffer,
+ const Array<Element> & elements,
+ SynchronizationTag tag) const;
+
+ virtual inline void unpackElementData(CommunicationBuffer & buffer,
+ const Array<Element> & elements,
+ SynchronizationTag tag);
+
+
+private:
+ /* ------------------------------------------------------------------------ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------ */
+ /// limit at which a point is considered as complitely broken
+ Real damage_limit;
+
+ /// internal pointer to the current damage vector
+ const Array<Real> * selected_damage;
+};
+
+
+#if defined (AKANTU_INCLUDE_INLINE_IMPL)
+# include "remove_damaged_weight_function_inline_impl.cc"
+#endif
+
+__END_AKANTU__
+
+#endif /* __AKANTU_REMOVE_DAMAGED_WEIGHT_FUNCTION_HH__ */
diff --git a/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.cc b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.cc
new file mode 100644
index 000000000..12d700ddf
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.cc
@@ -0,0 +1,100 @@
+/**
+ * @file remove_damaged_weight_function_inline_impl.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief Implementation of inline function of remove damaged weight function
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline void RemoveDamagedWeightFunction<spatial_dimension>::selectType(__attribute__((unused)) ElementType type1,
+ __attribute__((unused)) GhostType ghost_type1,
+ ElementType type2,
+ GhostType ghost_type2) {
+ /// select the damage array for a given type: For optimization
+ selected_damage = &(this->material.template getArray<Real>("damage", type2, ghost_type2));
+ // selected_damage = &mat.getDamage(type2, ghost_type2);
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline Real RemoveDamagedWeightFunction<spatial_dimension>::operator()(Real r,
+ const __attribute__((unused)) QuadraturePoint & q1,
+ const QuadraturePoint & q2) {
+ /// compute the weight
+ UInt quad = q2.global_num;
+
+ if(q1 == q2) return 1.;
+
+ Real D = (*selected_damage)(quad);
+ Real w = 0.;
+ if(D < damage_limit) {
+ Real alpha = std::max(0., 1. - r*r / this->R2);
+ w = alpha * alpha;
+ }
+ return w;
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline UInt RemoveDamagedWeightFunction<spatial_dimension>::getNbDataForElements(const Array<Element> & elements,
+ SynchronizationTag tag) const {
+ if(tag == _gst_mnl_weight)
+ return this->material.getModel().getNbQuadraturePoints(elements) * sizeof(Real);
+
+ return 0;
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline void RemoveDamagedWeightFunction<spatial_dimension>::packElementData(CommunicationBuffer & buffer,
+ const Array<Element> & elements,
+ SynchronizationTag tag) const {
+ if(tag == _gst_mnl_weight) {
+ ElementTypeMapArray<Real> & damage = this->material.template getInternal<Real>("damage");
+ this->material.packElementDataHelper(damage,
+ buffer,
+ elements);
+ }
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline void RemoveDamagedWeightFunction<spatial_dimension>::unpackElementData(CommunicationBuffer & buffer,
+ const Array<Element> & elements,
+ SynchronizationTag tag) {
+ if(tag == _gst_mnl_weight) {
+ ElementTypeMapArray<Real> & damage = this->material.template getInternal<Real>("damage");
+ this->material.unpackElementDataHelper(damage,
+ buffer,
+ elements);
+ }
+}
+
+
+
+
diff --git a/src/model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function.hh b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function.hh
new file mode 100644
index 000000000..8df6280e4
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function.hh
@@ -0,0 +1,85 @@
+/**
+ * @file remove_damaged_with damage_rate_weight_function.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief Removed damaged weight function for non local materials
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * Akantu is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "base_weight_function.hh"
+/* -------------------------------------------------------------------------- */
+#ifndef __AKANTU_REMOVE_DAMAGED_WITH_DAMAGE_RATE_WEIGHT_FUNCTION_HH__
+#define __AKANTU_REMOVE_DAMAGED_WITH_DAMAGE_RATE_WEIGHT_FUNCTION_HH__
+
+__BEGIN_AKANTU__
+/* -------------------------------------------------------------------------- */
+/* Remove damaged with damage rate weight function */
+/* -------------------------------------------------------------------------- */
+
+template<UInt spatial_dimension>
+class RemoveDamagedWithDamageRateWeightFunction : public BaseWeightFunction<spatial_dimension> {
+public:
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+ RemoveDamagedWithDamageRateWeightFunction(Material & material) : BaseWeightFunction<spatial_dimension>(material, "remove_damage_with_damage_rate") {
+ this->registerParam("damage_limit", this->damage_limit_with_damage_rate, 1, _pat_parsable, "Damage Threshold");
+ }
+
+ /* -------------------------------------------------------------------------- */
+ /* Base Weight Function inherited methods */
+ /* -------------------------------------------------------------------------- */
+ inline void selectType(__attribute__((unused)) ElementType type1,
+ __attribute__((unused)) GhostType ghost_type1,
+ ElementType type2,
+ GhostType ghost_type2);
+
+ inline Real operator()(Real r,
+ const __attribute__((unused)) QuadraturePoint & q1,
+ const QuadraturePoint & q2);
+
+ /* ------------------------------------------------------------------------ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------ */
+private:
+ /// limit at which a point is considered as complitely broken
+ Real damage_limit_with_damage_rate;
+
+ /// internal pointer to the current damage vector
+ const Array<Real> * selected_damage_with_damage_rate;
+
+ /// internal pointer to the current damage rate vector
+ const Array<Real> * selected_damage_rate_with_damage_rate;
+};
+
+#if defined (AKANTU_INCLUDE_INLINE_IMPL)
+# include "remove_damaged_with_damage_rate_weight_function_inline_impl.cc"
+#endif
+
+__END_AKANTU__
+
+#endif /* __AKANTU_REMOVE_DAMAGED_WITH_DAMAGE_WEIGHT_FUNCTION_HH__ */
diff --git a/src/model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function_inline_impl.cc b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function_inline_impl.cc
new file mode 100644
index 000000000..6f0eb805b
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function_inline_impl.cc
@@ -0,0 +1,65 @@
+/**
+ * @file remove_damaged_with damage_rate_weight_function_inline_impl.cc
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief Implementation of inline function of remove damaged with
+ * damage rate weight function
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline void RemoveDamagedWithDamageRateWeightFunction<spatial_dimension>::selectType(__attribute__((unused)) ElementType type1,
+ __attribute__((unused)) GhostType ghost_type1,
+ ElementType type2,
+ GhostType ghost_type2) {
+ /// select the damage arrays for given types: For optimization
+ selected_damage_with_damage_rate = &(this->material.template getArray<Real>("damage",type2, ghost_type2));
+ selected_damage_rate_with_damage_rate = &(this->material.template getArray<Real>("damage-rate",type2, ghost_type2));
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+inline Real RemoveDamagedWithDamageRateWeightFunction<spatial_dimension>::operator()(Real r,
+ const __attribute__((unused)) QuadraturePoint & q1,
+ const QuadraturePoint & q2) {
+ /// compute the weight
+ UInt quad = q2.global_num;
+
+ if(q1.global_num == quad) return 1.;
+
+ Real D = (*selected_damage_with_damage_rate)(quad);
+ Real w = 0.;
+ Real alphaexp = 1.;
+ Real betaexp = 2.;
+ if(D < damage_limit_with_damage_rate) {
+ Real alpha = std::max(0., 1. - pow((r*r / this->R2),alphaexp));
+ w = pow(alpha, betaexp);
+ }
+
+ return w;
+}
+
diff --git a/src/model/solid_mechanics/materials/weight_functions/stress_based_weight_function.hh b/src/model/solid_mechanics/materials/weight_functions/stress_based_weight_function.hh
new file mode 100644
index 000000000..e846802de
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/stress_based_weight_function.hh
@@ -0,0 +1,115 @@
+/**
+ * @file stress_based_weight_function.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief Removed damaged weight function for non local materials
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * Akantu is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "base_weight_function.hh"
+/* -------------------------------------------------------------------------- */
+#ifndef __AKANTU_STRESS_BASED_WEIGHT_FUNCTION_HH__
+#define __AKANTU_STRESS_BASED_WEIGHT_FUNCTION_HH__
+
+__BEGIN_AKANTU__
+/* -------------------------------------------------------------------------- */
+/* Stress Based Weight */
+/* -------------------------------------------------------------------------- */
+/// based on based on Giry et al.: Stress-based nonlocal damage model,
+/// IJSS, 48, 2011
+template<UInt spatial_dimension>
+class StressBasedWeightFunction : public BaseWeightFunction<spatial_dimension> {
+public:
+ /* ------------------------------------------------------------------------ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------ */
+ StressBasedWeightFunction(Material & material);
+
+ /* -------------------------------------------------------------------------- */
+ /* Base Weight Function inherited methods */
+ /* -------------------------------------------------------------------------- */
+ void init();
+
+ virtual inline void updateInternals(__attribute__((unused)) const ElementTypeMapArray<Real> & quadrature_points_coordinates);
+
+ void updatePrincipalStress(GhostType ghost_type);
+
+ inline void updateQuadraturePointsCoordinates(ElementTypeMapArray<Real> & quadrature_points_coordinates);
+
+ inline void selectType(ElementType type1, GhostType ghost_type1,
+ ElementType type2, GhostType ghost_type2);
+
+ inline Real operator()(Real r,
+ const QuadraturePoint & q1,
+ const QuadraturePoint & q2);
+
+ /// computation of ellipsoid
+ inline Real computeRhoSquare(Real r,
+ Vector<Real> & eigs,
+ Matrix<Real> & eigenvects,
+ Vector<Real> & x_s);
+
+private:
+
+ /* ------------------------------------------------------------------------ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------ */
+ /// tensile strength
+ Real ft;
+
+ /// prinicipal stresses
+ InternalField<Real> stress_diag;
+
+ /// for preselection of types (optimization)
+ Array<Real> * selected_stress_diag;
+
+ /// principal directions
+ InternalField<Real> stress_base;
+
+ /// for preselection of types (optimization)
+ Array<Real> * selected_stress_base;
+
+ // InternalField<Real> quadrature_points_coordinates;
+ /// for preselection of types (optimization)
+ Array<Real> * selected_position_1;
+ Array<Real> * selected_position_2;
+
+ /// lenght intrinisic to the material
+ InternalField<Real> characteristic_size;
+
+ /// for preselection of types (optimization)
+ Array<Real> * selected_characteristic_size;
+};
+
+#include "stress_based_weight_function_tmpl.hh"
+#if defined (AKANTU_INCLUDE_INLINE_IMPL)
+# include "stress_based_weight_function_inline_impl.cc"
+#endif
+
+__END_AKANTU__
+
+#endif /* __AKANTU_STRESS_BASED_WEIGHT_FUNCTION_HH__ */
diff --git a/src/model/solid_mechanics/materials/weight_function_tmpl.hh b/src/model/solid_mechanics/materials/weight_functions/stress_based_weight_function_inline_impl.cc
similarity index 61%
rename from src/model/solid_mechanics/materials/weight_function_tmpl.hh
rename to src/model/solid_mechanics/materials/weight_functions/stress_based_weight_function_inline_impl.cc
index 92dc9dabd..469892c30 100644
--- a/src/model/solid_mechanics/materials/weight_function_tmpl.hh
+++ b/src/model/solid_mechanics/materials/weight_functions/stress_based_weight_function_inline_impl.cc
@@ -1,279 +1,203 @@
/**
- * @file weight_function_tmpl.hh
+ * @file stress_based_weight_function_inline_impl.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
+ * @author Cyprien Wolff <cyprien.wolff@epfl.ch>
*
* @date creation: Fri Apr 13 2012
* @date last modification: Thu Jun 05 2014
*
- * @brief implementation of the weight function classes
+ * @brief Implementation of inline function of remove damaged with
+ * damage rate weight function
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
-
-/* -------------------------------------------------------------------------- */
-/* Stress based weight function */
-/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
-StressBasedWeightFunction<spatial_dimension>::StressBasedWeightFunction(Material & material) :
- BaseWeightFunction<spatial_dimension>(material, "stress_based"),
- stress_diag("stress_diag", material), selected_stress_diag(NULL),
- stress_base("stress_base", material), selected_stress_base(NULL),
- characteristic_size("lc", material), selected_characteristic_size(NULL) {
-
- this->registerParam("ft", this->ft, 0., _pat_parsable, "Tensile strength");
- stress_diag.initialize(spatial_dimension);
- stress_base.initialize(spatial_dimension * spatial_dimension);
- characteristic_size.initialize(1);
-}
-
-/* -------------------------------------------------------------------------- */
-template<UInt spatial_dimension>
-void StressBasedWeightFunction<spatial_dimension>::init() {
- const Mesh & mesh = this->material.getModel().getFEEngine().getMesh();
- for (UInt g = _not_ghost; g <= _ghost; ++g) {
- GhostType gt = GhostType(g);
- Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt);
- Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, gt);
- for(; it != last_type; ++it) {
- UInt nb_quadrature_points =
- this->material.getModel().getFEEngine().getNbQuadraturePoints(*it, gt);
- const Array<UInt> & element_filter = this->material.getElementFilter(*it, gt);
- UInt nb_element = element_filter.getSize();
-
- Array<Real> ones(nb_element*nb_quadrature_points, 1, 1.);
- Array<Real> & lc = characteristic_size(*it, gt);
- this->material.getModel().getFEEngine().integrateOnQuadraturePoints(ones,
- lc,
- 1,
- *it,
- gt,
- element_filter);
-
- for (UInt q = 0; q < nb_quadrature_points * nb_element; q++) {
- lc(q) = pow(lc(q), 1./ Real(spatial_dimension));
- }
- }
- }
-}
-
-
-/* -------------------------------------------------------------------------- */
-template<UInt spatial_dimension>
-void StressBasedWeightFunction<spatial_dimension>::updatePrincipalStress(GhostType ghost_type) {
- AKANTU_DEBUG_IN();
-
- const Mesh & mesh = this->material.getModel().getFEEngine().getMesh();
-
- Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
- Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
- for(; it != last_type; ++it) {
- Array<Real>::const_matrix_iterator sigma =
- this->material.getStress(*it, ghost_type).begin(spatial_dimension, spatial_dimension);
- Array<Real>::vector_iterator eigenvalues =
- stress_diag(*it, ghost_type).begin(spatial_dimension);
- Array<Real>::vector_iterator eigenvalues_end =
- stress_diag(*it, ghost_type).end(spatial_dimension);
- Array<Real>::matrix_iterator eigenvector =
- stress_base(*it, ghost_type).begin(spatial_dimension, spatial_dimension);
-
-#ifndef __trick__
- Array<Real>::iterator<Real> cl = characteristic_size(*it, ghost_type).begin();
-#endif
- UInt q = 0;
- for(;eigenvalues != eigenvalues_end; ++sigma, ++eigenvalues, ++eigenvector, ++cl, ++q) {
- sigma->eig(*eigenvalues, *eigenvector);
- *eigenvalues /= ft;
-#ifndef __trick__
- // specify a lower bound for principal stress based on the size of the element
- for (UInt i = 0; i < spatial_dimension; ++i) {
- (*eigenvalues)(i) = std::max(*cl / this->R, (*eigenvalues)(i));
- }
-#endif
- }
- }
- AKANTU_DEBUG_OUT();
+inline void StressBasedWeightFunction<spatial_dimension>::updateInternals(__attribute__((unused)) const ElementTypeMapArray<Real> & quadrature_points_coordinates) {
+ updatePrincipalStress(_not_ghost);
+ updatePrincipalStress(_ghost);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
inline void StressBasedWeightFunction<spatial_dimension>::selectType(ElementType type1,
GhostType ghost_type1,
ElementType type2,
GhostType ghost_type2) {
selected_stress_diag = &stress_diag(type2, ghost_type2);
selected_stress_base = &stress_base(type2, ghost_type2);
selected_characteristic_size = &characteristic_size(type1, ghost_type1);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
inline Real StressBasedWeightFunction<spatial_dimension>::operator()(Real r,
const QuadraturePoint & q1,
const QuadraturePoint & q2) {
Real zero = std::numeric_limits<Real>::epsilon();
if(r < zero) return 1.; // means x and s are the same points
const Vector<Real> & x = q1.getPosition();
const Vector<Real> & s = q2.getPosition();
Vector<Real> eigs =
selected_stress_diag->begin(spatial_dimension)[q2.global_num];
Matrix<Real> eigenvects =
selected_stress_base->begin(spatial_dimension, spatial_dimension)[q2.global_num];
Real min_rho_lc = selected_characteristic_size->begin()[q1.global_num];
Vector<Real> x_s(spatial_dimension);
x_s = x;
x_s -= s;
Real rho_2 = computeRhoSquare(r, eigs, eigenvects, x_s);
Real rho_lc_2 = std::max(this->R2 * rho_2, min_rho_lc*min_rho_lc);
// Real w = std::max(0., 1. - r*r / rho_lc_2);
// w = w*w;
Real w = exp(- 2*2*r*r / rho_lc_2);
return w;
}
/* -------------------------------------------------------------------------- */
template<>
inline Real StressBasedWeightFunction<1>::computeRhoSquare(__attribute__ ((unused)) Real r,
Vector<Real> & eigs,
__attribute__ ((unused)) Matrix<Real> & eigenvects,
__attribute__ ((unused)) Vector<Real> & x_s) {
return eigs[0];
}
/* -------------------------------------------------------------------------- */
template<>
inline Real StressBasedWeightFunction<2>::computeRhoSquare(__attribute__ ((unused)) Real r,
Vector<Real> & eigs,
Matrix<Real> & eigenvects,
Vector<Real> & x_s) {
Vector<Real> u1(eigenvects.storage(), 2);
Real cos_t = x_s.dot(u1) / (x_s.norm() * u1.norm());
Real cos_t_2;
Real sin_t_2;
Real sigma1_2 = eigs[0]*eigs[0];
Real sigma2_2 = eigs[1]*eigs[1];
#ifdef __trick__
Real zero = std::numeric_limits<Real>::epsilon();
if(std::abs(cos_t) < zero) {
cos_t_2 = 0;
sin_t_2 = 1;
} else {
cos_t_2 = cos_t * cos_t;
sin_t_2 = (1 - cos_t_2);
}
Real rhop1 = std::max(0., cos_t_2 / sigma1_2);
Real rhop2 = std::max(0., sin_t_2 / sigma2_2);
#else
cos_t_2 = cos_t * cos_t;
sin_t_2 = (1 - cos_t_2);
Real rhop1 = cos_t_2 / sigma1_2;
Real rhop2 = sin_t_2 / sigma2_2;
#endif
return 1./ (rhop1 + rhop2);
}
/* -------------------------------------------------------------------------- */
template<>
inline Real StressBasedWeightFunction<3>::computeRhoSquare(Real r,
Vector<Real> & eigs,
Matrix<Real> & eigenvects,
Vector<Real> & x_s) {
Vector<Real> u1(eigenvects.storage() + 0*3, 3);
//Vector<Real> u2(eigenvects.storage() + 1*3, 3);
Vector<Real> u3(eigenvects.storage() + 2*3, 3);
Real zero = std::numeric_limits<Real>::epsilon();
Vector<Real> tmp(3);
tmp.crossProduct(x_s, u3);
Vector<Real> u3_C_x_s_C_u3(3);
u3_C_x_s_C_u3.crossProduct(u3, tmp);
Real norm_u3_C_x_s_C_u3 = u3_C_x_s_C_u3.norm();
Real cos_t = 0.;
if(std::abs(norm_u3_C_x_s_C_u3) > zero) {
Real inv_norm_u3_C_x_s_C_u3 = 1. / norm_u3_C_x_s_C_u3;
cos_t = u1.dot(u3_C_x_s_C_u3) * inv_norm_u3_C_x_s_C_u3;
}
Real cos_p = u3.dot(x_s) / r;
Real cos_t_2;
Real sin_t_2;
Real cos_p_2;
Real sin_p_2;
Real sigma1_2 = eigs[0]*eigs[0];
Real sigma2_2 = eigs[1]*eigs[1];
Real sigma3_2 = eigs[2]*eigs[2];
#ifdef __trick__
if(std::abs(cos_t) < zero) {
cos_t_2 = 0;
sin_t_2 = 1;
} else {
cos_t_2 = cos_t * cos_t;
sin_t_2 = (1 - cos_t_2);
}
if(std::abs(cos_p) < zero) {
cos_p_2 = 0;
sin_p_2 = 1;
} else {
cos_p_2 = cos_p * cos_p;
sin_p_2 = (1 - cos_p_2);
}
Real rhop1 = std::max(0., sin_p_2 * cos_t_2 / sigma1_2);
Real rhop2 = std::max(0., sin_p_2 * sin_t_2 / sigma2_2);
Real rhop3 = std::max(0., cos_p_2 / sigma3_2);
#else
cos_t_2 = cos_t * cos_t;
sin_t_2 = (1 - cos_t_2);
cos_p_2 = cos_p * cos_p;
sin_p_2 = (1 - cos_p_2);
Real rhop1 = sin_p_2 * cos_t_2 / sigma1_2;
Real rhop2 = sin_p_2 * sin_t_2 / sigma2_2;
Real rhop3 = cos_p_2 / sigma3_2;
#endif
return 1./ (rhop1 + rhop2 + rhop3);
}
+
diff --git a/src/model/solid_mechanics/materials/weight_functions/stress_based_weight_function_tmpl.hh b/src/model/solid_mechanics/materials/weight_functions/stress_based_weight_function_tmpl.hh
new file mode 100644
index 000000000..30392c858
--- /dev/null
+++ b/src/model/solid_mechanics/materials/weight_functions/stress_based_weight_function_tmpl.hh
@@ -0,0 +1,116 @@
+/**
+ * @file stress_based_weight_function_tmpl.hh
+ *
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
+ *
+ * @date creation: Fri Apr 13 2012
+ * @date last modification: Thu Jun 05 2014
+ *
+ * @brief implementation of the stres based weight function classes
+ *
+ * @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 <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+StressBasedWeightFunction<spatial_dimension>::StressBasedWeightFunction(Material & material) :
+ BaseWeightFunction<spatial_dimension>(material, "stress_based"),
+ stress_diag("stress_diag", material), selected_stress_diag(NULL),
+ stress_base("stress_base", material), selected_stress_base(NULL),
+ characteristic_size("lc", material), selected_characteristic_size(NULL) {
+
+ this->registerParam("ft", this->ft, 0., _pat_parsable, "Tensile strength");
+ stress_diag.initialize(spatial_dimension);
+ stress_base.initialize(spatial_dimension * spatial_dimension);
+ characteristic_size.initialize(1);
+}
+
+/* -------------------------------------------------------------------------- */
+/// During intialization the characteristic sizes for all quadrature
+/// points are computed
+template<UInt spatial_dimension>
+void StressBasedWeightFunction<spatial_dimension>::init() {
+ const Mesh & mesh = this->material.getModel().getFEEngine().getMesh();
+ for (UInt g = _not_ghost; g <= _ghost; ++g) {
+ GhostType gt = GhostType(g);
+ Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt);
+ Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, gt);
+ for(; it != last_type; ++it) {
+ UInt nb_quadrature_points =
+ this->material.getModel().getFEEngine().getNbQuadraturePoints(*it, gt);
+ const Array<UInt> & element_filter = this->material.getElementFilter(*it, gt);
+ UInt nb_element = element_filter.getSize();
+
+ Array<Real> ones(nb_element*nb_quadrature_points, 1, 1.);
+ Array<Real> & lc = characteristic_size(*it, gt);
+ this->material.getModel().getFEEngine().integrateOnQuadraturePoints(ones,
+ lc,
+ 1,
+ *it,
+ gt,
+ element_filter);
+
+ for (UInt q = 0; q < nb_quadrature_points * nb_element; q++) {
+ lc(q) = pow(lc(q), 1./ Real(spatial_dimension));
+ }
+ }
+ }
+}
+
+
+/* -------------------------------------------------------------------------- */
+/// computation of principals stresses and principal directions
+template<UInt spatial_dimension>
+void StressBasedWeightFunction<spatial_dimension>::updatePrincipalStress(GhostType ghost_type) {
+ AKANTU_DEBUG_IN();
+
+ const Mesh & mesh = this->material.getModel().getFEEngine().getMesh();
+
+ Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
+ Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
+ for(; it != last_type; ++it) {
+ Array<Real>::const_matrix_iterator sigma =
+ this->material.getStress(*it, ghost_type).begin(spatial_dimension, spatial_dimension);
+ Array<Real>::vector_iterator eigenvalues =
+ stress_diag(*it, ghost_type).begin(spatial_dimension);
+ Array<Real>::vector_iterator eigenvalues_end =
+ stress_diag(*it, ghost_type).end(spatial_dimension);
+ Array<Real>::matrix_iterator eigenvector =
+ stress_base(*it, ghost_type).begin(spatial_dimension, spatial_dimension);
+
+#ifndef __trick__
+ Array<Real>::iterator<Real> cl = characteristic_size(*it, ghost_type).begin();
+#endif
+ UInt q = 0;
+ for(;eigenvalues != eigenvalues_end; ++sigma, ++eigenvalues, ++eigenvector, ++cl, ++q) {
+ sigma->eig(*eigenvalues, *eigenvector);
+ *eigenvalues /= ft;
+#ifndef __trick__
+ // specify a lower bound for principal stress based on the size of the element
+ for (UInt i = 0; i < spatial_dimension; ++i) {
+ (*eigenvalues)(i) = std::max(*cl / this->R, (*eigenvalues)(i));
+ }
+#endif
+ }
+ }
+ AKANTU_DEBUG_OUT();
+}
+
+