Commit 0b56d25e authored by Kevin Moreau's avatar Kevin Moreau

last version of code + cleanup



git-svn-id: https://svn.ec-nantes.fr/eXlibris/Applis/DamageBandDyn@1762 fbbead7c-fb4d-4173-aa67-51132c73c120
parent 3c35fff0
......@@ -121,8 +121,8 @@ public:
std::map<AOMD::mVertex*, AOMD::mEntity*>::iterator endv=vertex_map.end();
penetration_map.clear();
for(; itv!=endv; ++itv) {
mVertex* n = itv->first;
mEntity* e = itv->second;
AOMD::mVertex* n = itv->first;
AOMD::mEntity* e = itv->second;
xfem::xGeomElem geom_elem(e);
geom_elem.setUVWForXYZ(n->point());
std::pair<double,xfem::xVector> res;
......@@ -136,7 +136,7 @@ public:
for (; itp!=endp; ++itp) {
double penetration=itp->second.first;
if (penetration<0.000001) {
mVertex* n = itp->first;
AOMD::mVertex* n = itp->first;
xfem::xVector normal=itp->second.second;
std::vector<xfem::xValue<double>*> disp_vals;
std::vector<xfem::xValue<double>*> velo_vals;
......@@ -183,7 +183,7 @@ public:
for (; itp!=endp; ++itp) {
double penetration=itp->second.first;
if (penetration<0.000001) {
mVertex* n = itp->first;
AOMD::mVertex* n = itp->first;
xfem::xVector normal=itp->second.second;
std::vector<xfem::xValue<double>*> accel_vals;
xfem::xFiniteElement FEM;
......
......@@ -30,6 +30,18 @@ void PostProcessing::exportOnSpace(const std::string export_name, const int step
}
}
void PostProcessing::exportOnSpace(const std::string export_name, const int step, const xEval<double>& eval, const xIntegrationRule& integration_rule,
xIter begin,
xIter end)
{
if (export_manager.toExport(export_name, step, ""))
{
pexport.openFile(export_manager.getFileName());
Export(eval, pexport, export_name, integration_rule, begin, end);
pexport.closeFile();
}
}
void PostProcessing::exportOnSpace(const std::string export_name, const int step, const xEval<double>& eval, const xIntegrationRule& integration_rule,
boost::filter_iterator<xEntityFilter, xIter> begin,
boost::filter_iterator<xEntityFilter, xIter> end)
......
......@@ -33,6 +33,9 @@ public:
void reinit(xfem::xMesh* mesh);
void exportOnSpace(const std::string export_name, const int step, const xfem::xLevelSet& ls);
void exportOnSpace(const std::string export_name, const int step, const xfem::xEval<double>& eval, const xfem::xIntegrationRule& integration_rule,
xfem::xIter begin,
xfem::xIter end);
void exportOnSpace(const std::string export_name, const int step, const xfem::xEval<double>& eval, const xfem::xIntegrationRule& integration_rule,
boost::filter_iterator<xfem::xEntityFilter, xfem::xIter> begin,
boost::filter_iterator<xfem::xEntityFilter, xfem::xIter> end);
......
......@@ -14,6 +14,7 @@
#include "xVectorLevelSet.h"
#endif
using namespace AOMD;
using namespace xfem;
class LocalGeom
......
This diff is collapsed.
......@@ -26,6 +26,43 @@ class oOctree;
class oRefinementCriteria;
}
namespace oct_xfem_interf
{
const int init_physical_index = 100;
struct iClassifyCriteria
{
iClassifyCriteria() {}
virtual ~iClassifyCriteria() {}
virtual int operator()(oct::oOctree::const_iterator cell, int level, const int * ijk, oct::oOctree::const_iterator children_beg, oct::oOctree::const_iterator children_end) const = 0;
};
void InterfaceOctreeToAOMD(const oct::oKeyManager& key_manager, const oct::oOctree& octree,
xfem::xMesh& mesh,
const int id_offset,
boost::function<double (Trellis_Util::mPoint)> ptdt,
boost::function<double (Trellis_Util::mPoint)> ptdn,
std::map<int,int>& duplicate_nodes,
const bool simplex = true,
iClassifyCriteria* ptr_clfCriteria = NULL);
void InterfaceOctreeToAOMDNoHanging(const oct::oKeyManager& key_manager, const oct::oOctree& octree,
xfem::xMesh& mesh,
const bool coarser = false, const bool simplex = true, iClassifyCriteria* ptr_clfCriteria = NULL);
void OctreeToAOMDProjection(const oct::oKeyManager& key_manager, const oct::oField& field,
const std::map<int,int>& duplicate_nodes, const std::vector<double>& ls_values,
xfem::xMesh& mesh, xfem::xLevelSet& ls, const int id_offset);
void OctreeToAOMDNoHangingProjection(const oct::oKeyManager& key_manager, const oct::oField& field,
xfem::xMesh& mesh, xfem::xLevelSet& ls);
void InterfaceAOMDToOctree(xfem::xMesh& mesh, const xfem::xLevelSet& ls,
const std::map<int,int>& duplicate_nodes, std::vector<double>& ls_values,
const oct::oKeyManager& key_manager, oct::oField& field, const int id_offset);
}
class InterfaceOctreeAOMD
{
public:
......
......@@ -97,7 +97,6 @@ reInitLevelSetFromFrontFiltered::reInitLevelSetFromFrontFiltered(const std::vect
void reInitLevelSetFromFrontFiltered::visit(xLevelSet& ls, xRegion target)
{
#ifdef HAVE_CGAL
#ifndef OPTIMISATIONREINITLEVELSET
typedef xFilteredRegion<xIter, xEntityFilter>::FilterIter FilterIter;
double d=0.;
xPhysSurfByTagging domain(ls);
......@@ -123,29 +122,6 @@ void reInitLevelSetFromFrontFiltered::visit(xLevelSet& ls, xRegion target)
}
}
}
#else
typedef xFilteredRegion<xIter, xEntityFilter>::FilterIter FilterIter;
std::ostringstream oss;
oss << "tag_nodes_reInitLevelSetFromFrontFiltered_"<<this;
const unsigned int tag_nodes = AOMD_Util::Instance()->lookupMeshDataId(oss.str().c_str());
double d=0.;
for (std::vector<std::pair<xEntityFilter, xEntityFilter> >::const_iterator itf=filters.begin(); itf!=filters.end(); ++itf)
{
xFilteredRegion<xIter, xEntityFilter> fr_bnd(target.begin(2), target.end(2), itf->first);
xFilteredRegion<xIter, xEntityFilter> fr(target.begin(0), target.end(0), itf->second);
geom::xLevelSetTriToCGALEdgeAndTag generator(ls, tag_nodes);
geom::xDistanceNearestPoint<FilterIter, geom::xLevelSetTriToCGALEdgeAndTag> distance_nearest_point(fr_bnd.begin(), fr_bnd.end(), generator);
if (generator.isListNonEmpty())
{
for (FilterIter it=fr.begin(); it!=fr.end(); ++it)
{
d = distance_nearest_point.distance(static_cast<mVertex*>(*it));
if (ls(*it)>=0.) ls(*it)=d;
else ls(*it)=-d;
}
}
}
#endif
#endif
return;
}
......
......@@ -18,7 +18,9 @@
#include "tools.h"
unsigned int SpatialPartitioning::critical_frequency_tag = AOMD_Util::Instance()->lookupMeshDataId("critical_frequency_tag");
using namespace AOMD;
unsigned int SpatialPartitioning::critical_frequency_tag = AOMD_Util::Instance()->lookupMeshDataId("critical_frequency_tag");
unsigned int SpatialPartitioning::correction_frequency_tag = AOMD_Util::Instance()->lookupMeshDataId("correction_frequency_tag");
unsigned int SpatialPartitioning::resolution_frequency_tag = AOMD_Util::Instance()->lookupMeshDataId("resolution_frequency_tag");
unsigned int SpatialPartitioning::prediction_frequency_tag = AOMD_Util::Instance()->lookupMeshDataId("prediction_frequency_tag");
......
......@@ -187,19 +187,19 @@ private:
{
return (int)pow(2.,level-1);
}
inline void getElementNodes(mEntity* element, std::vector<mEntity*> & vector)
inline void getElementNodes(AOMD::mEntity* element, std::vector<AOMD::mEntity*> & vector)
{
for (int i = 0; i < element->size(0); ++i)
{
mEntity* node = element->get(0, i);
AOMD::mEntity* node = element->get(0, i);
vector.push_back(node);
}
}
inline void getNodeElements(mEntity* node, std::vector<mEntity*> & vector)
inline void getNodeElements(AOMD::mEntity* node, std::vector<AOMD::mEntity*> & vector)
{
for (int i = 0; i < node->size(dim); ++i)
{
mEntity* element = node->get(dim,i);
AOMD::mEntity* element = node->get(dim,i);
if (f_domain(element)) vector.push_back(element);
}
}
......
This diff is collapsed.
......@@ -134,13 +134,9 @@ private:
Observer& observer;
PostProcessing& post_processing;
std::set<AOMD::mEntity*> narrow_band;
std::set<AOMD::mEntity*> damage_band;
std::map<AOMD::mEntity*, AOMD::mEntity*> narrow_band_nodes;
std::map<AOMD::mEntity*, AOMD::mEntity*> damage_band_nodes;
const unsigned int narrow_tag;
const unsigned int damage_tag;
xfem::xSubMesh* narrow_band;
unsigned int narrow_band_attach_tag;
unsigned int damage_band_attach_tag;
};
#endif
......@@ -64,7 +64,6 @@
// XoctreeInterface
#if defined HAVE_OCTREE
#include "AdaptOctreeToAOMD.h"
#include "InterfaceOctreeToAOMDkev.h"
#include "oExport.h"
#endif
// xTLS
......
......@@ -25,15 +25,14 @@
#define WEIGHTEDVERSION 1
#define OPTIMISATIONPOINTWISE
// #define OPTIMISATIONREINITLEVELSET
// define: dead zone is part of computation
#define DEADZONE
// #define DEADZONE
// define: Takes into account values of the level set
// on boundaries. Coded for TLSSolverMode
#ifndef MODEBND
#define MODEBND
// #define MODEBND
#endif
// define: Use vector level set. Coded for TLSSolverMode
......
This diff is collapsed.
#ifndef _FilterManager_h
#define _FilterManager_h
#include "xEntityFilter.h"
#include "xRegion.h"
namespace xfem
{
class xPhysSurfByTagging;
class xIntegrationRuleBasic;
class xIntegrationRulePartition;
}
using namespace xfem;
enum FilterEnum
{
DOMAINGEOSUPPORT, DOMAININSUPPORT, DOMAINOUTSUPPORT, DOMAINGEO, DOMAININ,
DOMAINOUT, SANEMATTER, DAMAGEBAND, DAMAGEBANDSUPPORT, UNDAMAGED, UNDAMAGEDSUPPORT,
//HALFDAMAGEBAND,
INISOLC,
LINEARRAMPHEAVISIDE, VALUEONRAMPHEAVISIDE , DOMAINOUTSUPPORTBOUNDARY,
DOMAININSUPPORTELEMENTWISE, ALL
};
enum IntegrationRuleEnum
{
MASS, STIFFNESS, NEUMANN, EXPORT, CAUCHY, CONTACT
};
/*
* A filter manager, when domain are updated the information is
* udpated naturally on filters
*/
class FilterManager
{
public:
FilterManager(const xPhysSurfByTagging& domain_geo_);
FilterManager(const xPhysSurfByTagging& domain_geo_, const xPhysSurfByTagging& domain_in_, const xPhysSurfByTagging& domain_out_);
const xEntityFilter& getFilter(const FilterEnum filter_enum_);
const string getFilterName(const FilterEnum filter_enum_);
private:
xEntityFilter f_domain_geo;// omega
xEntityFilter f_domain_in;// gamma_0
xEntityFilter f_domain_out;// gamma_lc
xEntityFilter f_sane_matter;// omega \ (d=1)
xEntityFilter f_damage_band;// omega_0c
xEntityFilter f_damage_band_support;
xEntityFilter f_undamaged;
xEntityFilter f_undamaged_support;
//xEntityFilter f_half_damage_band;
xEntityFilter f_in_iso_lc;
xEntityFilter f_linear_ramp_heaviside;
xEntityFilter f_value_on_ramp_heaviside;
xEntityFilter f_domain_out_support_boundary;
xEntityFilter f_domain_in_support_element_wise;
xEntityFilter f_all;
};
/*
* A xFilteredRegion manager.
*/
class FilteredRegionManager
{
public:
typedef boost::filter_iterator<xEntityFilter, xIter> FilterIter;
typedef std::map<const FilterEnum, FilterIter*> Container;
typedef std::map<const FilterEnum, FilterIter*>::const_iterator ContainerIter;
FilteredRegionManager(const xRegion& all_, FilterManager& filter_manager_);
FilterIter* begin(const FilterEnum filter_enum_, const int nodal_=1);
FilterIter* end(const FilterEnum filter_enum_, const int nodal_=1);
void clear();
void exportGmsh(const FilterEnum filter_enum_, const int nodal_=1);
private:
xRegion all;
FilterManager& filter_manager;
Container container_begin;
Container container_end;
Container container_begin_0;
Container container_end_0;
};
/*
* An xIntegrationRule manager
*/
class IntegrationRuleManager
{
public:
typedef std::pair<const FilterEnum, const IntegrationRuleEnum> Pair;
typedef std::map<const Pair, xIntegrationRulePartition*> PartitionContainer;
typedef std::map<const Pair, xIntegrationRulePartition*>::const_iterator PartitionContainerIter;
typedef std::map<const IntegrationRuleEnum, xIntegrationRuleBasic*> BasicContainer;
typedef std::map<const IntegrationRuleEnum, xIntegrationRuleBasic*>::const_iterator BasicContainerIter;
IntegrationRuleManager(FilterManager& filter_manager_);
xIntegrationRulePartition* getIntegrationRule(const FilterEnum filter_enum_, const IntegrationRuleEnum integration_rule_enum_);
xIntegrationRuleBasic* getIntegrationRule(const IntegrationRuleEnum integration_rule_enum_);
void exportGmsh(const FilterEnum filter_enum_, const xRegion& all_);
private:
int getDegree(const IntegrationRuleEnum integration_rule_enum_);
FilterManager& filter_manager;
int mass_degree;
int stiffness_degree;
int neumann_degree;
int cauchy_degree;
int export_degree;
int contact_degree;
PartitionContainer partition_container;
BasicContainer basic_container;
};
#endif
This diff is collapsed.
#ifndef _Material_h
#define _Material_h
#include "xMaterial.h"
using namespace xfem;
class PBX : virtual public xMaterial
{
public:
PBX();
void checkProperties();
void sensitivityTo(const std::string& phys_token, xfem::xTensor4& sensitivity);
void sensitivityTo(const std::string& phys_token, xfem::xTensor4Isotropic& sensitivity);
void sensitivityTo(const std::string& phys_token, double& sensitivity);
void sensitivityTo(const std::string& phys_token, std::string& sensitivity);
void computeCurrentState();
static inline double delta(int i, int j) {return ((i == j)?1.0:0.0);}
static void SetElasticStiffnessIsotropic (double E, double nu, xfem::xTensor4& tensor);
static void SetElasticComplianceIsotropic(double E, double nu, xTensor4& compliance);
protected:
double lam, mu;
double E, nu;
xfem::xTensor4 elastic_stiffness;
xfem::xTensor4 elastic_compliance;
};
class CompressibleNeoHookean : public xMaterial
{
public:
CompressibleNeoHookean();
void sensitivityTo(const std::string& phys_token, double& sensitivity);
void sensitivityTo(const std::string& phys_token, xfem::xTensor4& sensitivity);
void computeCurrentState();
double computeStrainNRJDensity();
void checkProperties();
inline double delta(int i, int j)
{
return ((i == j)?1.0:0.0);
}
private:
double lambda, mu;
xTensor2 ident2;
};
class UnidimensionalElastic : public xMaterial
{
public:
UnidimensionalElastic();
void checkProperties();
void sensitivityTo(const std::string& phys_token, double& sensitivity);
private:
};
#endif
#ifndef _UnidimAlgorithm_h
#define _UnidimAlgorithm_h
#include "xAssembler.h"
#include "xData.h"
#include "xField.h"
#include "xIntegrationRule.h"
#include "Export.h"
#include "InputManager.h"
#include "UnidimCommandOnGeomElem.h"
using namespace xfem;
/*
* Hacked version of TreatmentOfNatEnv for unidim
*/
void TreatmentOfNatEnv1D(const xField& listFunctionSpace, xAssembler& assembler, const xIntegrationRulePartition& integrator, xData* data, xBoundary& groups, const double& t, const xEntityFilter filter=xAcceptAll())
{
xParseData& main_info = InputManagerSingleton::instance();
double intensity_left = main_info.getDouble("LOAD_INTENSITY_LEFT");
double intensity_right = main_info.getDouble("LOAD_INTENSITY_RIGHT");
double rate = main_info.getDouble("LOAD_RATE");
for (xPhysicalEnv::iterator it = data->PhysicalEnv->begin(); it != data->PhysicalEnv->end(); ++it) {
const xEnv& env = *it;
if (env.Phys == "TRACTION_X") {
assert(env.Type == FIX);
double val = 0.;
if (env.Phys == "TRACTION_X") {
if (it==data->PhysicalEnv->begin()) {
val = intensity_left;
} else {
val = intensity_right;
}
}
map<double,double> evolution_map;
if (rate!=0) evolution_map.insert(pair<double,double>(0.,0.));
evolution_map.insert(pair<double,double>(rate,1.));
evolution_map.insert(pair<double,double>(100.,1.));
xPieceWiseLinearDoubleToDouble evolution(evolution_map);
xEvalConstant<double> constant(val);
EvalSpaceFunctionCrossTimeFunction<double> traction(constant,evolution);
traction.setTime(t);
xFormLinearWithLoad<xValOperator<xIdentity<double> >, EvalSpaceFunctionCrossTimeFunction<double> > lin(traction);
xClassRegion bc(data->mesh, env.Entity, env.getDimension());
if (it==data->PhysicalEnv->begin()) {
xExportSensors& sensor_manager = ExportSensorSingleton::instance();
sensor_manager.measureAll(traction, t, "TreatmentOfNatEnv1D_loading");
}
xFilteredRegion<xClassIter, xEntityFilter> fr_bc(bc.begin(), bc.end(), filter);
Assemble(lin, assembler, integrator, listFunctionSpace, fr_bc.begin(), fr_bc.end(), xUpperAdjacency());
}
}
}
/*
* Shift level sets
*/
inline void UpdateLevelSets(const double delta_length, xLevelSet& ls_in_, xLevelSet& ls_out_)
{
xShiftLevelSetModifier shift_ls_modifier(delta_length);
ls_in_.accept(shift_ls_modifier);
ls_out_.accept(shift_ls_modifier);
}
inline void ExportForMatplotlib(const vector<double>& x, const vector<double>& y, const char* file_name)
{
std::fstream filestr(file_name, std::fstream::out);
std::vector<double>::const_iterator it_x=x.begin();
std::vector<double>::const_iterator it_y=y.begin();
for (;it_x!=x.end();++it_x,++it_y) {
filestr << *it_x << " " << *it_y << std::endl;
}
filestr.close();
}
template <class ITER>
void Export1D(const xEval<double>& eval, const char* file_name, const xIntegrationRule& integration_rule, ITER it, ITER end, bool on_integration_point=false)
{
if (on_integration_point) {
MatplotlibIntegrationPoint plot(eval, file_name);
ApplyCommandOnIntegrationRule(plot, integration_rule, it, end);
} else {
MatplotlibPoint plot(eval, file_name);
ApplyCommandOnIntegrationRule(plot, integration_rule, it, end);
}
}
inline void Export1D(const char* file_name, const xLevelSet& ls)
{
std::fstream filestr(file_name, std::fstream::out);
if (filestr.is_open()) {
xRegion region = ls.getSupport();
for (xIter it=region.begin(0); it!=region.end(0); ++it) {
mVertex* v = (mVertex*) *it;
mPoint p = v->point();
filestr << p(0) << " " << ls(v) << std::endl;
}
}
filestr.close();
}
/*
* Rebuild domains
*/
inline void UpdateDomain(xPhysSurfByTagging& domain_geo, xPhysSurfByTagging& domain_in, xPhysSurfByTagging& domain_out, FilteredRegionManager& fr_manager)
{
domain_geo.clean();
domain_in.clean();
domain_out.clean();
domain_geo.update(false, false, false);
domain_in.update(false, true, false);
domain_out.update(false, true, false);
fr_manager.clear();
}
#endif
#ifndef _UnidimCommandOnGeomElem_h
#define _UnidimCommandOnGeomElem_h
#include "xCommandOnGeomElem.h"
#include "xEval.h"
using namespace xfem;
/*
* Export dof values to .txt file
*/
class MatplotlibPoint : public xCommandOnGeomElem {
public:
MatplotlibPoint(const xEval<double>& eval_, const char* file_name) :
eval(eval_),
filestr(file_name, fstream::out)
{
}
virtual ~MatplotlibPoint()
{
filestr.close();
}
void execute(xGeomElem* geom_integ)
{
mEntity* elt = geom_integ->getEntity();
for (int i=0; i<elt->size(0); ++i) {
geom_integ->setUVWForVertex(i);
mVertex* v = (mVertex*) elt->get(0,i);
int sub = 0;
if (geom_appro->getEntity()!=geom_integ->getEntity()) {
geom_appro->setUVWForXYZ(geom_integ->getXYZ());
sub = 1;
}
else geom_appro->setUVW(geom_integ->getUVW());
mPoint p = geom_integ->getXYZ();
double val;
eval(geom_appro, geom_integ, val);
filestr << p(0) << " " << val << " " << sub << " " << v->getId() << endl;
}
}
private:
const xEval<double>& eval;
fstream filestr;
};
/*
* Export integration point values to .txt file
*/
class MatplotlibIntegrationPoint : public xCommandOnGeomElem {
public:
MatplotlibIntegrationPoint(const xEval<double>& eval_, const char* file_name) :
eval(eval_),
filestr(file_name, fstream::out)
{
}
virtual ~MatplotlibIntegrationPoint()
{
filestr.close();
}
void execute(xGeomElem* geom_integ)
{
for (int i=0; i<geom_integ->GetNbIntegrationPoints(); ++i) {
geom_integ->setUVW(i);
if (geom_appro->getEntity()!=geom_integ->getEntity()) geom_appro->setUVWForXYZ(geom_integ->getXYZ());
else geom_appro->setUVW(geom_integ->getUVW());
mPoint p_integ = geom_integ->getXYZ();
double val;
eval(geom_appro, geom_integ, val);
filestr << p_integ(0) << " " << val << endl;
}
}
private:
const xEval<double>& eval;
fstream filestr;
};
/*
*
*/
// class PlasticRadialReturn : public xCommandOnGeomElem {
// public:
// PlasticRadialReturn(const xEval<double>& eval_strain_, xFieldPointwise& isotropic_hardening_pointwise_field_, xFieldPointwise& isotropic_hardening_pointwise_old_field_) : eval_strain(eval_strain_), isotropic_hardening_pointwise_field(isotropic_hardening_pointwise_field_), isotropic_hardening_pointwise_old_field(isotropic_hardening_pointwise_old_field_) {}
// ~PlasticRadialReturn() {}
// void execute(xGeomElem* geom_integ) {
// for (int i=0; i<geom_integ->GetNbIntegrationPoints(); ++i) {
// geom_integ->setUVW(i);
// if (geom_appro->getEntity()!=geom_integ->getEntity()) {
// geom_appro->setUVWForXYZ(geom_integ->getXYZ());
// }
// else geom_appro->setUVW(geom_integ->getUVW());
// // mPoint p1 = geom_appro->getXYZ();
// // mPoint p2 = geom_integ->getXYZ();
// // cout << p1(0) << " " << p2(0) << endl;
// double plastic_threshold_val = std::numeric_limits<double>::max();
// double plastic_threshold_val_0 = 0.;
// double plastic_threshold_val_derivative = 0.;
// double delta_lambda = 0.;
// double tolerance = 1e-4;
// xEvalField<xIdentity<double> ,xFieldPointwise> eval_plastic_cumulate_strain(isotropic_hardening_pointwise_field);
// int counter = 0;
// while (counter < 10)
// {
// xEvalConstant<double> eval_delta_lambda(delta_lambda);
// xEvalBinary<std::plus<double> > eval_plastic_cumulate_strain_newton(eval_plastic_cumulate_strain, eval_delta_lambda);
// EvalPlasticStrain eval_plastic_strain(eval_plastic_cumulate_strain_newton);
// EvalElasticStrain eval_elastic_strain(eval_strain, eval_plastic_strain);
// EvalSurface eval_surface;
// EvalEffectiveStress eval_effective_stress(eval_elastic_strain, eval_surface);
// EvalIsotropicHardening eval_isotropic_hardening(eval_plastic_cumulate_strain_newton, eval_surface);
// EvalIsotropicHardeningDerivative eval_isotropic_hardening_derivative(eval_plastic_cumulate_strain, eval_surface);
// EvalPlasticThresholdFunction eval_plastic_threshold_function(eval_effective_stress, eval_isotropic_hardening);
// EvalPlasticThresholdFunctionDerivative eval_plastic_threshold_function_derivative(eval_isotropic_hardening_derivative, eval_surface);
// eval_plastic_threshold_function(geom_appro, geom_integ, plastic_threshold_val);
// if (plastic_threshold_val==0.) plastic_threshold_val_0 = plastic_threshold_val;
// if (plastic_threshold_val<=tolerance*fabs(plastic_threshold_val_0)) break;
// eval_plastic_threshold_function_derivative(geom_appro, geom_integ, plastic_threshold_val_derivative);
// delta_lambda += - plastic_threshold_val/plastic_threshold_val_derivative;
// counter++;
// message(__func__,"plastic threshold val", plastic_threshold_val);
// message(__func__,"iteration radial return "+toString(counter));
// message(__func__,"delta lambda", delta_lambda);
// }
// double plastic_cumulate_strain = 0.;
// eval_plastic_cumulate_strain(geom_appro, geom_integ, plastic_cumulate_strain);
// isotropic_hardening_pointwise_old_field.setVal(geom_appro, geom_integ, plastic_cumulate_strain);