Commit 90d823ac authored by Kévin Moreau's avatar Kévin Moreau

This commit bring all this changes

* add an epsilon_ratio (between 0 and 1) in order to control the fmeik internal epsilon involved when comparing
scalar values (a < b + epsilon is now a <b + epsilon_ratio*"edge length")

it affects:
  - FastMarchingInterface.h
  - TLSSolver.{h,cc}
  - FormulationQSRemeshAniso.{h,cc}

* modify anticipated end of program: now lambda doesn't have to belong to [0, 1].

it affects:
  - FormulationQS.cc

* change dissipated energy calculation. It is based on the primitive of function h(d).

it affects:
  - FormulationQS.cc
  - MaterialElasticDamage.cc
  - MaterialFunction.{h,cc}

* compute nonlocal part bounding box

it affects:
  - TLSGeom.{h,cc}
  - TLSSolver.cc

* minor simplification when transfering phi field, now only "int_duplicated" region is used.

it affects:
  - MeshGeneration.{h,cc}
  - TLSGeom.cc
  - Formulation.cc
  - TLSSolver.cc

* bug correction when transfer boundary and interior "in" after remeshing.

it affects:
  - TLSGeom.cc

* bug correction in correctLevelSetField: add a forgotten absolute value.

it affects:
  - TLSSolver.cc

* place guards with error message to help user detect bugs

it affects:
  - FormulationQS.cc
  - TLSSolver.cc
parent 4ee7604c
......@@ -5,11 +5,7 @@
#ifndef _FastMarchingInterface_h_
#define _FastMarchingInterface_h_
// Xfem
#include "xField.h"
#include "xMesh.h"
#include "xRegion.h"
// Xinterfaces
#include "xDenseMatrix.h"
//
#ifdef HAVE_FASTMARCHING
#include "meshinterfacexRegion.h"
......@@ -59,6 +55,7 @@ void FastMarchingReinit(std::function<double (AOMD::mVertex*)> get_val,
ITERTRIALKNOWN trial_known_begin, ITERTRIALKNOWN trial_known_end,
ITERTOFIND tofind_begin, ITERTOFIND tofind_end,
ITERKNOWN known_it,
double epsilon_ratio,
std::function<double (const AOMD::mVertex&)> f_func=[](const AOMD::mVertex&){return 1.;}) {
#ifndef HAVE_FASTMARCHING
std::cout << "FastMarchingReinit requires HAVE_FASTMARCHING defined to work" << std::endl;
......@@ -77,7 +74,7 @@ void FastMarchingReinit(std::function<double (AOMD::mVertex*)> get_val,
entitytovertexiteratorconvertor<ITERTRIALKNOWN> trial_known_conv_begin(trial_known_begin);
entitytovertexiteratorconvertor<ITERTRIALKNOWN> trial_known_conv_end(trial_known_end);
entitytovertexiteratorconvertor<ITERKNOWN> known_conv_it(known_it);
fmeik(mi, ls, trial_known_conv_begin, trial_known_conv_end, known_conv_it, f_func, gls);
fmeik(mi, ls, trial_known_conv_begin, trial_known_conv_end, known_conv_it, f_func, epsilon_ratio, gls);
for(auto it=region.begin(0); it!=region.end(0); ++it) {
auto v=static_cast<AOMD::mVertex*>(*it);
double val;
......@@ -88,13 +85,15 @@ void FastMarchingReinit(std::function<double (AOMD::mVertex*)> get_val,
}
template<typename ITERKNOWN,
typename ITERTOFIND>
typename ITERTOFIND,
typename MATRIX>
void FastMarchingModeExtension(std::function<double (AOMD::mVertex*)> get_val,
std::function<void (AOMD::mVertex*, double)> set_val,
const xfem::xRegion& region,
ITERKNOWN known_begin, ITERKNOWN known_end,
ITERTOFIND tofind_begin, ITERTOFIND tofind_end,
lalg::xDenseMatrix& coeffs,
MATRIX& coeffs,
double epsilon_ratio,
const std::function<double (const AOMD::mVertex&)>& f_func=[](const AOMD::mVertex&){return 1.;}) {
#ifndef HAVE_FASTMARCHING
std::cout << "FastMarchingModeExtension requires HAVE_FASTMARCHING defined to work" << std::endl;
......@@ -121,7 +120,7 @@ void FastMarchingModeExtension(std::function<double (AOMD::mVertex*)> get_val,
std::vector<const AOMD::mVertex*> vec;
vec.reserve(nb_modes);
entitytovertexiteratorconvertor<std::back_insert_iterator<std::vector<const AOMD::mVertex*>>> dummy_it(std::back_inserter(vec));
fmeik(mi, ls, known_conv_begin, known_conv_end, dummy_it, f_func, gls, vn);
fmeik(mi, ls, known_conv_begin, known_conv_end, dummy_it, f_func, epsilon_ratio, gls, vn);
int j=0;
for(auto it=region.begin(0); it!=region.end(0); ++it, ++j) {
auto v=static_cast<AOMD::mVertex*>(*it);
......
......@@ -426,19 +426,19 @@ bool QSFormulation::exportStep(int step) {
tls_solver.exportStep(step);
if(load_factor*ref_disp.mag()>0.1*geom.getBodyCharacteristicLength()) {
std::cout<<"Warning: anticipated end of program since maximum displacement is greater than 10\% of body characteristic length."<<std::endl;
std::cout<<"It happens when body is cut into two or more pieces!"<<std::endl;
return true;
}
else if(load_factor>1.) {
std::cout<<"Warning: anticipated end of program since computation of load factor is wrong."<<std::endl;
return true;
}
else if(std::fabs(load_factor_variation)>1.) {
std::cout<<"Warning: anticipated end of program since computation of load factor variation is wrong."<<std::endl;
return true;
}
// if(load_factor*ref_disp.mag()>0.1*geom.getBodyCharacteristicLength()) {
// std::cout<<"Warning: anticipated end of program since maximum displacement is greater than 10\% of body characteristic length."<<std::endl;
// std::cout<<"It happens when body is cut into two or more pieces!"<<std::endl;
// return true;
// }
// else if(load_factor>1.) {
// std::cout<<"Warning: anticipated end of program since computation of load factor is wrong."<<std::endl;
// return true;
// }
// else if(std::fabs(load_factor_variation)>1.) {
// std::cout<<"Warning: anticipated end of program since computation of load factor variation is wrong."<<std::endl;
// return true;
// }
NonUniformMaterialSensitivity_c<double> eval_delta_dissipated_energy("delta_dissipated_energy", variab_manager);
double delta_dissipated_energy=0.;
......@@ -446,7 +446,13 @@ bool QSFormulation::exportStep(int step) {
ApplyCommandOnIntegrationRule(integ_command, geom.getIntegRuleSmart(), begin, end);
post_pro.exportOnTime("delta_dissipated_energy", (double)step, delta_dissipated_energy);
dissipated_energy+=delta_dissipated_energy;
// dissipated_energy+=delta_dissipated_energy;
NonUniformMaterialSensitivity_c<double> eval_dissipated_energy("dissipated_energy", variab_manager);
dissipated_energy=0.;
xIntegrateEvalCommand<xEval<double> > integ_command2(eval_dissipated_energy, dissipated_energy);
ApplyCommandOnIntegrationRule(integ_command2, geom.getIntegRuleSmart(), begin, end);
post_pro.exportOnTime("dissipated_energy", (double)step, dissipated_energy);
post_pro.exportOnTime("load_factor", (double)step, load_factor);
......
......@@ -55,6 +55,7 @@ QSFormulRemeshAniso::QSFormulRemeshAniso(TLSGeom& geom, TLSSolver& tls_solver,
const xData& data, const xParseData& parse_data,
PreProcessing& pre_pro, PostProcessing& post_pro, std::function<void (xMesh&)> tagger) :
QSFormulation(geom, tls_solver, data, parse_data, pre_pro, post_pro),
epsilon_ratio(static_cast<bool>(parse_data.getInt("do_bnd_in"))?1.e-6:0.),
tagger(tagger){}
......@@ -291,14 +292,14 @@ void QSFormulRemeshAniso::buildLSAdaptationField(){
targetRegion,
sm_bnd_nonlocal.begin(0), sm_bnd_nonlocal.end(0),
sm_local.begin(0), sm_local.end(0),
std::back_inserter(dummy));
std::back_inserter(dummy), epsilon_ratio);
FastMarchingReinit([](mVertex* v){ return 0.; },
[this](mVertex* v, double val){ this->lsetAdapt(v)=val; },
targetRegion,
sm_bnd_nonlocal.begin(0), sm_bnd_nonlocal.end(0),
sm_int_nonlocal.begin(0), sm_int_nonlocal.end(0),
std::back_inserter(dummy));
std::back_inserter(dummy), epsilon_ratio);
/* KEV: refactor using interface of TLSDuctile
meshinterfacexRegion mi(targetRegion);
......
......@@ -38,6 +38,7 @@ private:
void projectPhiOnMesh(xMesh *mesh_new, std::map<mEntity *, double> &projectedValues);
void xAssignPhiValues(xField &field, xMesh *mesh, std::map<mEntity *, double> &projectedValues);
std::list<AOMD::mEntity*> constructDelocalizedList(xfem::xLevelSet&);
const double epsilon_ratio;
};
......
......@@ -86,6 +86,9 @@ void ElastoDam::sensitivityTo(const std::string& phys_token, double& sensitivity
double delta_d=d-old->scalar("d");
sensitivity=Y_c*hardening_function->getVal(d)*delta_d;
}
else if (phys_token == "dissipated_energy") {
sensitivity=Y_c*hardening_function->getPrimitive(curr->scalar("d"));
}
else if(phys_token == "Y") {
sensitivityTo("a", sensitivity);
}
......
......@@ -32,6 +32,10 @@ double AffinePowerHardeningFunction::getVal(double x) const {
double AffinePowerHardeningFunction::getGrad(double x) const {
return p*pow(1.+a*x, p-1)*a;
}
double AffinePowerHardeningFunction::getPrimitive(double x) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
PolynomialHardeningFunction::PolynomialHardeningFunction(const std::string& filename) {
signature.register_table("COEFFICIENTS");
......@@ -62,6 +66,10 @@ double PolynomialHardeningFunction::getGrad(double x) const {
}
return p;
}
double PolynomialHardeningFunction::getPrimitive(double x) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
TabulatedHardeningFunction::TabulatedHardeningFunction(const std::string& filename) {
signature.register_string("VALUES_FILENAME");
......@@ -106,6 +114,10 @@ double TabulatedHardeningFunction::getGrad(double x) const {
y1=get<2>(f[i]);
return y0+(y1-y0)/(x1-x0)*(x-x0) ;
}
double TabulatedHardeningFunction::getPrimitive(double x) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
ExponentialHardeningFunction::ExponentialHardeningFunction(const std::string& filename) {
signature.register_scalar("HARDENING_COEFF");
......@@ -119,6 +131,10 @@ double ExponentialHardeningFunction::getVal(double x) const {
double ExponentialHardeningFunction::getGrad(double x) const {
return c*exp(c*x);
}
double ExponentialHardeningFunction::getPrimitive(double x) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
RationalHardeningFunction::RationalHardeningFunction(const std::string& filename) {
signature.register_scalar("HARDENING_COEFF");
......@@ -134,6 +150,10 @@ double RationalHardeningFunction::getVal(double x) const {
double RationalHardeningFunction::getGrad(double x) const {
return c*p/pow(1.-x, p+1);
}
double RationalHardeningFunction::getPrimitive(double x) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
LogarithmHardeningFunction::LogarithmHardeningFunction(const std::string& filename) {
signature.register_scalar("HARDENING_COEFF");
......@@ -147,6 +167,10 @@ double LogarithmHardeningFunction::getVal(double x) const {
double LogarithmHardeningFunction::getGrad(double x) const {
return c/(1.-x);
}
double LogarithmHardeningFunction::getPrimitive(double x) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
CustomHardeningFunction::CustomHardeningFunction(const std::string& filename) {
signature.register_scalar("HARDENING_RATIO");
......@@ -164,6 +188,10 @@ double CustomHardeningFunction::getVal(double x) const {
double CustomHardeningFunction::getGrad(double x) const {
return (1.-r)*c*exp(c*x)+r*p/pow(1.-x, p+1);
}
double CustomHardeningFunction::getPrimitive(double x) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
CohesiveLinearHardeningFunction::CohesiveLinearHardeningFunction(const std::string& filename) {
signature.register_scalar("CZM_WF");
......@@ -180,40 +208,52 @@ CohesiveLinearHardeningFunction::CohesiveLinearHardeningFunction(const std::stri
double E = properties.scalar("YOUNG_MODULUS");
Yc0=0.5*sigmac*sigmac/E;
lda=2.*sigmac*lc/(E*wf);
assert(lda<0.5);
}
double CohesiveLinearHardeningFunction::getVal(double d) const {
d=min(d,0.99999);
// d=min(d,0.99999);
double I(0.),dI(0.),S(0.),dS(0.);
// double I(0.),dI(0.),S(0.),dS(0.);
I=(2.-d)/sqrt(1.-d)-2;
dI=0.5*d*pow(1.-d,-1.5);
// I=(2.-d)/sqrt(1.-d)-2;
// dI=0.5*d*pow(1.-d,-1.5);
S=pow(1+lda*I,-1);
dS=-lda*dI*pow(1+lda*I,-2);
// S=pow(1+lda*I,-1);
// dS=-lda*dI*pow(1+lda*I,-2);
return Yc0 * (S*S*pow(1.-d,-2)+2.*dS*S*d*pow(1.-d,-1));
// return S*S*pow(1.-d,-2)+2.*dS*S*d*pow(1.-d,-1);
double A=sqrt(1.-d)+lda*pow(1.-sqrt(1.-d), 2.);
double dA=(lda-0.5)/sqrt(1.-d)-lda;
return 1./(A*A)-2.*d/(A*A*A)*dA;
}
double CohesiveLinearHardeningFunction::getGrad(double d) const {
d=min(d,0.99999);
// d=min(d,0.99999);
double I(0.),dI(0.),ddI(0.),S(0.),dS(0.),ddS(0.);
// double I(0.),dI(0.),ddI(0.),S(0.),dS(0.),ddS(0.);
I=(2.-d)/sqrt(1.-d)-2;
dI=0.5*d*pow(1.-d,-1.5);
ddI=0.5*(1.+0.5*d)*pow(1.-d,-2.5);
// I=(2.-d)/sqrt(1.-d)-2;
// dI=0.5*d*pow(1.-d,-1.5);
// ddI=0.5*(1.+0.5*d)*pow(1.-d,-2.5);
S=pow(1+lda*I,-1);
dS=-lda*dI*pow(1+lda*I,-2);
ddS=-lda*(ddI*(1+lda*I)-2*lda*dI*dI)*pow(1+lda*I,-3);
// S=pow(1+lda*I,-1);
// dS=-lda*dI*pow(1+lda*I,-2);
// ddS=-lda*(ddI*(1+lda*I)-2*lda*dI*dI)*pow(1+lda*I,-3);
// // return 2.*S*S*pow(1.-d,-3)+4.*dS*S*pow(1.-d,-2)+2.*d*dS*dS*pow(1.-d,-1)+2.*ddS*S*pow(1.-d,-1);
// return 2.*S*S*pow(1.-d,-3)+4.*dS*S*pow(1.-d,-2)+2.*d*dS*dS*pow(1.-d,-1)+2.*d*ddS*S*pow(1.-d,-1);
double A=sqrt(1.-d)+lda*pow(1.-sqrt(1.-d), 2.);
double dA=(lda-0.5)/sqrt(1.-d)-lda;
double d2A=(lda-0.5)/(2.*pow(1.-d, 1.5));
return -4.*dA/(A*A*A)+6.*d*dA/(A*A*A*A)-2.*d*d2A/(A*A*A);
}
return Yc0*(2.*S*S*pow(1.-d,-3)+4.*dS*S*pow(1.-d,-2)+2.*d*dS*dS*pow(1.-d,-1)+2.*ddS*S*pow(1.-d,-1));
double CohesiveLinearHardeningFunction::getPrimitive(double d) const {
return d/pow(sqrt(1.-d)+lda*pow(1.-sqrt(1.-d), 2.), 2.);
}
CohesiveBilinearHardeningFunction::CohesiveBilinearHardeningFunction(const std::string& filename) {
......@@ -245,7 +285,6 @@ CohesiveBilinearHardeningFunction::CohesiveBilinearHardeningFunction(const std::
sigmak_dimless=sigmak/sigmac;
wk_dimless=w1_dimless*(1-sigmak_dimless);
Yc0=0.5*sigmac*sigmac/E;
lda=2.*sigmac*lc/(E*wf);
a=wk_dimless/(lda*sigmak_dimless);
phik_dimless=0.5*a*(sqrt(1.+4./a)-1.);
......@@ -268,7 +307,7 @@ double CohesiveBilinearHardeningFunction::getVal(double d) const {
S=pow(C,alpha)*pow(1+C*lda*I,-1);
dS=-pow(C,alpha+1.)*lda*dI*pow(1+C*lda*I,-2);
return Yc0 * (S*S*pow(1.-d,-2)+2.*dS*S*d*pow(1.-d,-1));
return S*S*pow(1.-d,-2)+2.*dS*S*d*pow(1.-d,-1);
}
double CohesiveBilinearHardeningFunction::getGrad(double d) const {
......@@ -289,8 +328,13 @@ double CohesiveBilinearHardeningFunction::getGrad(double d) const {
dS=-pow(C,alpha+1.)*lda*dI*pow(1+C*lda*I,-2);
ddS=-pow(C,alpha+1.)*lda*(ddI*(1+C*lda*I)-2*C*lda*dI*dI)*pow(1+C*lda*I,-3);
return Yc0*(2.*S*S*pow(1.-d,-3)+4.*dS*S*pow(1.-d,-2)+2.*d*dS*dS*pow(1.-d,-1)+2.*ddS*S*pow(1.-d,-1));
return 2.*S*S*pow(1.-d,-3)+4.*dS*S*pow(1.-d,-2)+2.*d*dS*dS*pow(1.-d,-1)+2.*ddS*S*pow(1.-d,-1);
}
double CohesiveBilinearHardeningFunction::getPrimitive(double d) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
LinearCouplingFunction::LinearCouplingFunction(const std::string& filename) {
signature.register_scalar("P_CRIT");
......@@ -304,6 +348,11 @@ double LinearCouplingFunction::getVal(double d) const {
double LinearCouplingFunction::getGrad(double d) const {
return p_c;
}
double LinearCouplingFunction::getPrimitive(double d) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
ExponentialCouplingFunction::ExponentialCouplingFunction(const std::string& filename) {
signature.register_scalar("P_CRIT");
......@@ -317,3 +366,7 @@ double ExponentialCouplingFunction::getVal(double d) const {
double ExponentialCouplingFunction::getGrad(double d) const {
return p_c/(1.-std::min(0.999, d));
}
double ExponentialCouplingFunction::getPrimitive(double d) const {
std::cout<<"TODO"<<std::endl;
return 0.;
}
......@@ -15,6 +15,7 @@ public:
MaterialFunction();
virtual double getVal(double) const = 0 ;
virtual double getGrad(double) const = 0 ;
virtual double getPrimitive(double) const = 0 ;
protected:
xfem::xTensors properties;
xfem::xTensorsSignature signature;
......@@ -25,6 +26,7 @@ public:
AffinePowerHardeningFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
double a, p;
};
......@@ -42,6 +44,7 @@ public:
PolynomialHardeningFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
xTable coeff;
int deg;
......@@ -52,6 +55,7 @@ public:
TabulatedHardeningFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
vector<tuple<double,double,double>> f ;
};
......@@ -61,6 +65,7 @@ public:
ExponentialHardeningFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
double c;
};
......@@ -70,6 +75,7 @@ public:
RationalHardeningFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
double c, p;
};
......@@ -79,6 +85,7 @@ public:
LogarithmHardeningFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
double c;
};
......@@ -88,23 +95,27 @@ public:
CustomHardeningFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
double r, c, p;
};
// Equivalence with CZM hardening
// WARNING :
// - The sigma_c coefficient that is given inside this class must be coherent with
// the prescribed Y_CRIT given in *.mat file (Y_c=1/2 sigma_c^2 / E).
// - CohesiveLinearHardening implements the small h(d) from Parilla-Gomez paper.
// - LC has to be given in file *.mat
// - YC_CRIT must be set to 1 (it is then recalculated from CZM parameters, and
// taken into account into H(d)
// - taken into account into h(d)
class CohesiveLinearHardeningFunction : public MaterialFunction {
public:
CohesiveLinearHardeningFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
double lc, wf, sigmac;
double Yc0, lda ;
double lda ;
};
class CohesiveBilinearHardeningFunction : public MaterialFunction {
......@@ -112,9 +123,10 @@ public:
CohesiveBilinearHardeningFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
double lc, wf, w1, sigmak, sigmac;
double Yc0, wk_dimless, w1_dimless, sigmak_dimless, lda, a, phik_dimless, dk ;
double wk_dimless, w1_dimless, sigmak_dimless, lda, a, phik_dimless, dk ;
};
class LinearCouplingFunction : public MaterialFunction {
......@@ -122,6 +134,7 @@ public:
LinearCouplingFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
double p_c;
};
......@@ -131,6 +144,7 @@ public:
ExponentialCouplingFunction(const std::string&);
virtual double getVal(double) const;
virtual double getGrad(double) const;
virtual double getPrimitive(double) const;
private:
double p_c;
};
......
......@@ -365,13 +365,11 @@ void NodeDispatcher::operator()(mEntity* np) {
}
NodeDuplicator::NodeDuplicator(xMesh& comp_mesh,
xSubMesh& sm_duplicated,
const unsigned int was_created_by_tag,
const unsigned int is_crack_tag,
const unsigned int partition_tag,
const unsigned int duplicated_node_tag) :
comp_mesh(comp_mesh),
sm_duplicated(sm_duplicated),
dim(comp_mesh.dim()),
was_created_by_tag(was_created_by_tag),
is_crack_tag(is_crack_tag),
......
......@@ -409,7 +409,6 @@ private:
class NodeDuplicator {
public:
NodeDuplicator(xfem::xMesh&,
xfem::xSubMesh&,
const unsigned int,
const unsigned int,
const unsigned int,
......@@ -421,13 +420,11 @@ private:
void duplicateNode(AOMD::mEntity* np, Iter begin_fi, Iter end_fi) {
const Trellis_Util::mPoint& pt=static_cast<const AOMD::mVertex*>(np)->point();
AOMD::mEntity* nd=comp_mesh.createVertex(comp_mesh.getNewVertexId(), pt(0), pt(1), pt(2), np->getClassification());
sm_duplicated.add(nd);
nd->attachEntity(was_created_by_tag, np); // only required to propagate classification
std::for_each(begin_fi, end_fi, PartitionDispatcher(duplicated_node_tag, [nd](AOMD::mEntity*){return nd;}));
}
xfem::xMesh& comp_mesh;
xfem::xSubMesh& sm_duplicated;
const unsigned int dim;
const unsigned int was_created_by_tag;
const unsigned int is_crack_tag;
......
......@@ -99,6 +99,30 @@ public:
private:
xMesh& mesh;
};
void ComputeBoundingBox(const xSubMesh& dom, mPoint& min, mPoint& max) {
auto it=dom.begin(0);
auto ite=dom.end(0);
if(it==ite) {
for(int j=0; j<3; ++j) {
min(j)=0.;
max(j)=0.;
}
return;
}
min=static_cast<mVertex*>(*it)->point();
max=min;
++it;
for(; it!=ite; ++it) {
const auto& p=static_cast<mVertex*>(*it)->point();
for(int j=0; j<3; ++j) {
if(min(j)>p(j)) min(j)=p(j);
if(max(j)<p(j)) max(j)=p(j);
}
}
}
IntegrationRuleSmart::IntegrationRuleSmart() {}
......@@ -340,6 +364,7 @@ void TLSGeom::buildTLS() { // TODO the purpose of this is to limit phi to a part
}
}
void TLSGeom::buildNonlocal(const std::list<mEntity*>& delocalized_list) {
auto& sm_nonlocal=getMesh().createSubMesh("nonlocal");
for(auto e: delocalized_list) {
......@@ -349,6 +374,19 @@ void TLSGeom::buildNonlocal(const std::list<mEntity*>& delocalized_list) {
}
}
void TLSGeom::computeNonlocalBoundingBox() const {
auto& sm_nonlocal=getMesh().getSubMesh("nonlocal");
mPoint min, max;
ComputeBoundingBox(sm_nonlocal, min, max);
int step=Observer::tell<int>("step");
post_pro.exportOnTime("nonlocal_bb_min_x", step, min(0));
post_pro.exportOnTime("nonlocal_bb_min_y", step, min(1));
post_pro.exportOnTime("nonlocal_bb_min_z", step, min(2));
post_pro.exportOnTime("nonlocal_bb_max_x", step, max(0));
post_pro.exportOnTime("nonlocal_bb_max_y", step, max(1));
post_pro.exportOnTime("nonlocal_bb_max_z", step, max(2));
}
void TLSGeom::updateNonlocal(const std::list<mEntity*>& delocalized_list) {
auto& sm_nonlocal=getMesh().getSubMesh("nonlocal");
for(auto e: delocalized_list) {
......@@ -491,7 +529,7 @@ void TLSGeom::transferBndIntInNonlocal() {
}
for(auto v: todump_int) {
sm_int_in_nonlocal.del(v);
}
}
xEvalConstant<double> eval_one(1.);
xIntegrationRuleBasic integ_rule(0);
......@@ -782,7 +820,7 @@ void TLSGeom::buildCompMesh() {
auto& sm_duplicated=comp_mesh.createSubMesh("duplicated");
std::for_each(sm_bnd_fds.begin(0), sm_bnd_fds.end(0), NodeDispatcher(comp_mesh, was_created_by_tag, is_crack_tag, partition_tag, duplicated_node_tag));
std::for_each(sm_int_fds.begin(0), sm_int_fds.end(0), NodeDuplicator(comp_mesh, sm_duplicated, was_created_by_tag, is_crack_tag, partition_tag, duplicated_node_tag));
std::for_each(sm_int_fds.begin(0), sm_int_fds.end(0), NodeDuplicator(comp_mesh, was_created_by_tag, is_crack_tag, partition_tag, duplicated_node_tag));
std::for_each(sm_fds.begin(), sm_fds.end(), FaceDuplicator(comp_mesh, sm_duplicated, was_created_by_tag, is_crack_tag, partition_tag, duplicated_tag, duplicated_node_tag));
std::for_each(sm_fds.begin(), sm_fds.end(), DelUpdateAdj(comp_mesh));
std::for_each(sm_int_fds.begin(1), sm_int_fds.end(1), DelUpdateAdj(comp_mesh));
......
......@@ -56,6 +56,7 @@ public:
double getElementSize() const;
double getBodyCharacteristicLength() const;
void computeNonlocalBoundingBox() const;
void updateDomains(std::list<AOMD::mEntity*>&);
void transferDomains();
......
......@@ -6,7 +6,6 @@
#include "FastMarchingInterface.h"
#include "Import.h"
#include "Observer.h"
#include "TLSAlgorithm.h"
#include "TLSGeom.h"
#include "TLSSolver.h"
#include "Util.h"
......@@ -20,10 +19,13 @@
#include "ValueOldAndCurrent.h"
#include "xCSRMatrix.h"//for xAssembler
#include "xCSRVector.h"
#include "xDenseMatrix.h"
#include "xGenericSparseMatrix.h"
#include "xGenericSparseMatrixTraitPolicy.h"
#include "TLSAlgorithm.h"
using namespace xfem;
using namespace xtls;
using namespace lalg;
......@@ -54,7 +56,8 @@ TLSSolver::TLSSolver(TLSGeom& geom, const xParseData& parse_data,
damageShapeFunctions::StringToDamageShape(parse_data.getString("damage_shape")),
xAcceptAll(), xAcceptAll()),
do_crack_cut(static_cast<bool>(parse_data.getInt("do_crack_cut"))),
is_newly_delocalized(false) {
is_newly_delocalized(false),
epsilon_ratio(static_cast<bool>(parse_data.getInt("do_bnd_in"))?1.e-6:0.) {
space_factory->setSpaceProductOrder(1);
field.insert(space_factory->getSpace("fem"));
......@@ -178,7 +181,7 @@ void TLSSolver::reinitLevelSetField() {
nonlocal,
bnd_nonlocal.begin(0), bnd_nonlocal.end(0),
int_nonlocal.begin(0), int_nonlocal.end(0),
std::back_inserter(bnd_in_nonlocal_vec));
std::back_inserter(bnd_in_nonlocal_vec), epsilon_ratio);
correctLevelSetField();
geom.buildBndIntInNonlocal(bnd_in_nonlocal_vec);
geom.deleteBndIntNonlocal();
......@@ -194,7 +197,7 @@ void TLSSolver::reinitLevelSetField() {
void TLSSolver::correctLevelSetField() {
const double lc=parse_data.getDouble("lc");
const double epsilon=geom.getElementSize()*parse_data.getDouble("fit_to_vertices_ratio");
const double epsilon=fabs(geom.getElementSize()*parse_data.getDouble("fit_to_vertices_ratio"));
double corr_less=lc;
double corr_more=lc;
......@@ -277,25 +280,27 @@ void TLSSolver::delocalizeLevelSetField() {
delocalized_list.clear();
is_newly_delocalized=true;
}
geom.computeNonlocalBoundingBox();
}
void TLSSolver::transferLevelSetField() {
geom.transferDomains();
declareField("duplicated");
declareField("int_duplicated");
deleteState();
post_pro.exportOnSpace("transfered_before_phi", Observer::tell<int>("step"), eval_phi, geom.getIntegRuleBasic(0), geom.begin("tls"), geom.end("tls"));
std::vector<mEntity*> dummy;
FastMarchingReinit([this](mVertex* v){ std::vector<double> vals; this->field.getVals(v, vals); return vals[0]; },
[this](mVertex* v, double val){ this->field.setVal(v, val); },
geom.getRegion("duplicated"),
geom.begin(0, "bnd_fds"), geom.end(0, "bnd_fds"),
geom.begin(0, "int_duplicated"), geom.end(0, "int_duplicated"),
std::back_inserter(dummy));
// TransferField(field, field,
// [](xValue<double>* v1, xValue<double>* v2) { v2->setVal(v1->getVal()); },
// geom.begin(0, "duplicated"), geom.end(0, "duplicated"));
// An update of the iso-l_c discretisation is necessary to do this.
// std::vector<mEntity*> dummy;
// FastMarchingReinit([this](mVertex* v){ std::vector<double> vals; this->field.getVals(v, vals); return vals[0]; },
// [this](mVertex* v, double val){ this->field.setVal(v, val); },
// geom.getRegion("duplicated"),
// geom.begin(0, "bnd_fds"), geom.end(0, "bnd_fds"),
// geom.begin(0, "int_duplicated"), geom.end(0, "int_duplicated"),
// std::back_inserter(dummy), epsilon_ratio);
TransferField(field, field,
[](xValue<double>* v1, xValue<double>* v2) { v2->setVal(v1->getVal()); },
geom.begin(0, "int_duplicated"), geom.end(0, "int_duplicated"));
deleteField("int_fds");
declareState();
......@@ -327,7 +332,7 @@ void TLSSolver::delocalizeMeanField() {
nonlocal,
bnd_in_nonlocal.begin(0), bnd_in_nonlocal.end(0),
int_in_nonlocal.begin(0), int_in_nonlocal.end(0),
coeffs);
coeffs, epsilon_ratio);
geom.deleteBndIntInNonlocal();
if(nonlocality_matrix) {
......
......@@ -82,8 +82,6 @@ public:
void saveStep(int);
void exportStep(int);
void exportMeanField(std::string, int);// heavy
protected:
void transferDomains();
private:
void restoreStep(PreProcessing&);
void correctLevelSetField();
......@@ -116,6 +114,7 @@ private:
const bool do_crack_cut;
bool is_newly_delocalized;
const double epsilon_ratio;
};
#endif
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment