Commit 5b561468 authored by Alexis SALZMAN's avatar Alexis SALZMAN

Merge branch 'master' into para

Conflicts:
	src/TLSSolver.cc
        resolved : just two header set on the same line
parents 36828c9e 3415cdbc
Pipeline #220 skipped
#include "FastMarchingInterface.h"
#include "TLSAlgorithm.h"
#include "xAlgorithm.h"
#include "xLevelSet.h"
#include "xMesh.h"
......@@ -25,19 +23,25 @@ namespace lalg {
using namespace xfem;
using namespace lalg;
using namespace xtls;
using Trellis_Util::mPoint;
using AOMD::mEntity;
using AOMD::mVertex;
int main(int argc, char* argv[]) {
if(argc!=3) {
std::cout<<"First argument: .msh file; Second argument: a boolean to say if bnd_in is to do."<<std::endl;
std::abort();
}
xMesh mesh(argv[1]);
// domain creation
xPlane p1(mPoint(0.1,0.1,0.), xVector(1.,0.,0.));
xPlane p2(mPoint(0.205,0.205,0.), xVector(-1.,0.,0.));
xLevelSet level_set(&mesh, xUnion(p1, p2));
xGrownHalfLine grown(mPoint(0.5, 0.5, 0.), xVector(1.,0.,0.), 0.26);
xLevelSet level_set(&mesh, xCompl(grown));
// xPlane p1(mPoint(0.1,0.1,0.), xVector(1.,0.,0.));
// xPlane p2(mPoint(0.205,0.205,0.), xVector(-1.,0.,0.));
// xLevelSet level_set(&mesh, xUnion(p1, p2));
xPhysSurfByTagging phys_surf(level_set);
xEntityFilter filter=std::bind1st(std::mem_fun(&xPhysSurfByTagging::strictOut), &phys_surf);
xSubMesh domain("domain", mesh);
......@@ -69,21 +73,22 @@ int main(int argc, char* argv[]) {
}
}
xRegion all(&mesh);
xRegion region(&domain);
std::cout<<"region "<<region.size(2)<<std::endl;
std::cout<<"bnd "<<bnd.size(0)<<std::endl;
std::cout<<"interior "<<interior.size(0)<<std::endl;
level_set.load(p1);
// level_set.load(p1);
// phi field
xDoubleManager double_manager;
xSpacePolynomialLagrange space("space", xSpace::SCALAR, 1);
xField field(&double_manager, space);
DeclareInterpolation(field, xValueCreator<xValueDouble>(), region.begin(), region.end());
DeclareState(field, xStateDofCreator<>(double_manager, "dofs"), region.begin(), region.end());
for(xIter it=region.begin(0); it!=region.end(0); ++it) {
DeclareInterpolation(field, xValueCreator<xValueDouble>(), all.begin(), all.end());
DeclareState(field, xStateDofCreator<>(double_manager, "dofs"), all.begin(), all.end());
for(xIter it=all.begin(0); it!=all.end(0); ++it) {
xFiniteElementKeysOnly keys;
keys.setKeys(*it, field.begin(), field.end());
std::vector<xValue<double>*> vals;
......@@ -95,7 +100,7 @@ int main(int argc, char* argv[]) {
xEvalField<xIdentity<double> > eval_field(field);
xIntegrationRuleBasic integ_rule(1);
xExportGmshAscii pexport;
Export(eval_field, pexport, "phi_given", integ_rule, region.begin(), region.end());
Export(eval_field, pexport, "phi_given", integ_rule, all.begin(), all.end());
std::function<double (const AOMD::mVertex&)> f_func=[](const AOMD::mVertex&){ return 1.; };
......@@ -106,13 +111,20 @@ int main(int argc, char* argv[]) {
region,
bnd.begin(0), bnd.end(0),
interior.begin(0), interior.end(0),
std::back_inserter(bnd_in_vec));
std::back_inserter(bnd_in_vec), 1.e-6);
Export(eval_field, pexport, "phi", integ_rule, region.begin(), region.end());
Export(eval_field, pexport, "phi", integ_rule, all.begin(), all.end());
xSubMesh bnd_in("bnd_in", mesh);
for(auto v: bnd_in_vec) {
bnd_in.add(const_cast<AOMD::mEntity*>(v));
if(atoi(argv[2])>0) {
for(auto v: bnd_in_vec) {
bnd_in.add(const_cast<AOMD::mEntity*>(v));
}
}
else {
for(auto it=bnd.begin(0); it!=bnd.end(0); ++it) {
bnd_in.add(*it);
}
}
xSubMesh interior_in("interior_in", mesh);
xAcceptInSubMesh is_in_bnd_in(bnd_in);
......@@ -133,12 +145,20 @@ int main(int argc, char* argv[]) {
region,
bnd_in.begin(0), bnd_in.end(0),
interior_in.begin(0), interior_in.end(0),
coeffs);
coeffs, 1.e-6);
std::cout<<"nb modes "<<bnd_in.size(0)<<std::endl;
// modes export
for(int i=0; i<bnd_in.size(0); ++i) {
for(xIter it=all.begin(0); it!=all.end(0); ++it) {
xFiniteElementKeysOnly keys;
keys.setKeys(*it, field.begin(), field.end());
std::vector<xValue<double>*> vals;
vals.reserve(keys.sizeKey());
double_manager.getValPtr(keys.beginKey(), keys.endKey(), vals);
vals[0]->setVal(0.);
}
int j=0;
for(xIter it=region.begin(0); it!=region.end(0); ++it, ++j) {
xFiniteElementKeysOnly keys;
......@@ -150,7 +170,7 @@ int main(int argc, char* argv[]) {
}
std::ostringstream oss;
oss<<i;
Export(eval_field, pexport, "mode_"+oss.str(), integ_rule, region.begin(), region.end());
Export(eval_field, pexport, "mode_"+oss.str(), integ_rule, all.begin(), all.end());
}
int j=0;
......@@ -167,7 +187,7 @@ int main(int argc, char* argv[]) {
}
vals[0]->setVal(sum);
}
Export(eval_field, pexport, "sum", integ_rule, region.begin(), region.end());
Export(eval_field, pexport, "sum", integ_rule, all.begin(), all.end());
return 0;
}
......@@ -83,7 +83,6 @@ PostProcessing::PostProcessing(const xParseData& parse_data_, xExport &pexport_,
filter(filter) ,
pexport(pexport_) {
export_manager.init(parse_data, "export_manager");
createCurvesPlotFiles(parse_data.getString("curve_plot_file_ext"));
}
void PostProcessing::init(xMesh& mesh) {
......@@ -189,71 +188,3 @@ void PostProcessing::archiveSave(int step) {
system("rm -f *_save.*");
}
}
void PostProcessing::createCurvesPlotFiles(std::string type) {
if (type=="matlab") {
std::ofstream file;
// Generation of a matlab file for the plot
file.open("plot_global_response.m", ios::out);
file << "clear all" << endl;
file << "close all" << endl;
file << " " << endl;
file << "% Table extraction" << endl;
file << "table = importdata('sensor_val_.txt', ' ', 1);" << endl;
file << "lis = [0 ; table.data(2:end,1)];" << endl;
file << "lisU = max(min(lisU, 1.e5), -1.e5);" << endl;
file << "table = importdata('sensor_val_lambda.txt', ' ', 1);" << endl;
file << " " << endl;
file << "% List of values" << endl;
file << "lisF = [0 ; table.data(2:end,1)];" << endl;
file << "lisFx = [0 ; table.data(2:end,2)];" << endl;
file << "lisFy = [0 ; table.data(2:end,3)];" << endl;
file << "lisFz = [0 ; table.data(2:end,4)];" << endl;
file << "lisF = max(min(lisF, 1.e5), -1.e5);" << endl;
file << " " << endl;
file << "% Plot" << endl;
file << "h1 = figure();" << endl;
file << "set(h1, 'position', [200 600 400 300]);" << endl;
file << "hold on" << endl;
file << "plot(lisU, lisF, 'b-+');" << endl;
file << "legend('norm of F');" << endl;
file << " " << endl;
file << "h2 = figure();" << endl;
file << "set(h2, 'position', [200 200 400 300]);" << endl;
file << "hold on" << endl;
file << "plot(lisU, lisFx, 'b-+');" << endl;
file << "plot(lisU, lisFy, 'r--+');" << endl;
file << "plot(lisU, lisFz, 'g-.+');" << endl;
file << "legend('Fx','Fy','Fz');" << endl;
file << "hold off" << endl;
file << " " << endl;
file << " " << endl;
file << " " << endl;
file << "% Table extraction" << endl;
file << "table_lambda = importdata('sensor_val_lambda.txt', ' ', 1);" << endl;
file << "table_force_x = importdata('sensor_val_force_x.txt', ' ', 1);" << endl;
file << "table_force_y = importdata('sensor_val_force_y.txt', ' ', 1);" << endl;
file << " " << endl;
file << "% List of values" << endl;
file << "lis_lambda = [0 ; table_lambda.data(2:end,2)];" << endl;
file << "lis_force_x = [0 ; table_force_x.data(2:end,2)];" << endl;
file << "lis_force_y = [0 ; table_force_y.data(2:end,2)];" << endl;
file << "lis_force_x = max(min(lis_force_x, 1.e5), -1.e5);" << endl;
file << "lis_force_y = max(min(lis_force_y, 1.e5), -1.e5);" << endl;
file << " " << endl;
file << "% Plot" << endl;
file << "h3 = figure()" << endl;
file << "set(h3, 'position', [700 200 400 300]);" << endl;
file << "hold on" << endl;
file << "plot(lis_lambda, lis_force_x, 'b-+');" << endl;
file << "plot(lis_lambda, lis_force_y, 'r--+');" << endl;
file << "legend('Fx','Fy');" << endl;
file << "hold off" << endl;
file.close();
}
}
......@@ -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);
......
......@@ -84,7 +84,7 @@ void Formulation::transferDispField(double coeff) {
CustomValueDoubleSelector::which_val=i;
TransferField(disp_field, disp_field,
[size](xValue<double>* v1, xValue<double>* v2) { v2->setVal(v1->getVal()); },
geom.begin(0, "duplicated"), geom.end(0, "duplicated"));// TODO degree one only!!!
geom.begin(0, "int_duplicated"), geom.end(0, "int_duplicated"));// TODO degree one only!!!
}
CustomValueDoubleSelector::which_val=0;
deleteDispField("int_fds");
......
......@@ -33,7 +33,7 @@ public:
const xfem::xData&, const xParseData&,
PostProcessing&);
virtual ~Formulation();
public:
// Simple accessors
int getNbDofs() const;
std::string getMaterialClass() const;
......
......@@ -52,20 +52,17 @@ QSFormulation::QSFormulation(TLSGeom& geom, TLSSolver& tls_solver,
declareDispField();
applyEssentialEnv(1.);
declareDofs();
if(pre_pro.doRestart()) {
restoreStep(pre_pro);
}
declareDispMeasField();
for(std::string str: {"a", "b", "c", "alpha", "beta"}) {
tls_solver.registerMeanVector("mean_"+str);
}
if(getMaterialClass()=="elastic_damage_nonlinear") {
is_nonlinear=true;
geom.buildLinearAndNonLinear(tls_solver.getPhiEval());
}
if(pre_pro.doRestart()) {
restoreStep(pre_pro);
}
std::cout<<"delta phi max "<<delta_phi_max<<std::endl;
}
......@@ -161,10 +158,10 @@ void QSFormulation::computeLoadFactorVariation() {
for(int i=0; i<a.size(); ++i) {
double f=a[i]*load_factor*load_factor+b[i]*load_factor+c[i];
double load_factor_variation_trial=(alpha[i]*delta_phi_max-f)/beta[i];
double delta_phi_trial=std::max(0., f+beta[i]*load_factor_variation_trial)/alpha[i];
if(delta_phi_trial>0.) {
load_factor_variation=std::min(load_factor_variation_trial, load_factor_variation);
}
// double delta_phi_trial=std::max(0., f+beta[i]*load_factor_variation_trial)/alpha[i];
// if(delta_phi_trial>0.) {
load_factor_variation=std::min(load_factor_variation_trial, load_factor_variation);
// }
}
std::cout<<"load factor variation "<<load_factor_variation<<std::endl;
Observer::listen("load_factor_variation", load_factor_variation);
......@@ -247,11 +244,28 @@ void QSFormulation::computeDeltaLevelSetField() {
const auto& alpha=tls_solver.getMeanVector("mean_alpha");
const auto& beta=tls_solver.getMeanVector("mean_beta");
auto& delta_phi=tls_solver.getMeanVector("delta_phi");
double min_delta_phi=delta_phi_max;
double max_delta_phi=0.;
int cnt=0;
for(int i=0; i<delta_phi.size(); ++i) {
const double f=a[i]*load_factor*load_factor+b[i]*load_factor+c[i];
delta_phi[i]=std::max(0., f+beta[i]*load_factor_variation)/alpha[i];
const double buf=std::max(0., f+beta[i]*load_factor_variation)/alpha[i];
delta_phi[i]=buf;
min_delta_phi=std::min(buf, min_delta_phi);
max_delta_phi=std::max(buf, max_delta_phi);
if(buf>0.) {
++cnt;
}
}
post_pro.exportOnTime("nb_active_modes", Observer::tell<int>("step"), (double)std::count_if(delta_phi.begin(), delta_phi.end(), [](double a){ return a>0.; }));
std::cout<<"delta_phi min modal value "<<min_delta_phi<<" delta_phi modal max value "<<max_delta_phi<<std::endl;
if(min_delta_phi<0. || std::fabs(max_delta_phi-delta_phi_max)>std::numeric_limits<double>::epsilon()) {
std::cout<<"Error: min of delta_phi modal values must be in [0, delta_phi_max] AND max close to delta_phi_max!"<<std::endl;
std::cout<<" if min_delta_phi is negative, this is due to a negative mean_alpha modal value,"<<std::endl;
std::cout<<" which is due to a negative hardening_function grad or some negative mode."<<std::endl;
std::cout<<" Check the positivity of your hardening getGrad function."<<std::endl;
std::abort();
}
post_pro.exportOnTime("nb_active_modes", Observer::tell<int>("step"), (double)cnt);
}
void QSFormulation::updateLevelSetField() {
......@@ -431,14 +445,14 @@ bool QSFormulation::exportStep(int step) {
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;
}
// 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 +460,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);
......@@ -489,23 +509,17 @@ void QSFormulation::restoreStep(PreProcessing& pre_pro) {
updateMaterialVariables(tls_solver.getDamageEval());
writeStrain("eps_ref");
tls_solver.swapOldAndCurrentLevelSetField();
updateMaterialVariables(tls_solver.getDamageEval());
SmoothMaterialVariable_c<double> eval("d", variab_manager);
post_pro.exportOnSpace("restored_d", pre_pro.getStep(), eval, geom.getIntegRuleBasic(0), geom.begin(), geom.end());
// TODO: the following lines recreate the link between integ_mesh and comp_mesh.
// It is far from optimal.
// The best solution is to export all relations between integ_mesh and comp_mesh
// and restablish it.
if(tls_solver.isCrackable()) {
geom.buildFullyDamagedSupport(tls_solver.getPhiMinusLcEval());
tls_solver.cutCrack();
tls_solver.transferLevelSetField();
transferDispField();
transferMaterialVariables();
updateMaterialVariables(tls_solver.getDamageEval(), "duplicated");
cutCrack();
}
tls_solver.delocalizeMeanField();
updateMaterialVariables(tls_solver.getDamageEval());
SmoothMaterialVariable_c<double> eval("d", variab_manager);
post_pro.exportOnSpace("restored_d", pre_pro.getStep(), eval, geom.getIntegRuleBasic(0), geom.begin(), geom.end());
}
......@@ -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;
};
......
......@@ -26,25 +26,30 @@ QSFormulationSuperimposed::QSFormulationSuperimposed(TLSGeom& geom, TLSSolver& t
QSFormulationSuperimposed::~QSFormulationSuperimposed() {}
void QSFormulationSuperimposed::applyEssentialEnv(double scale) {
xClassRegion bc(&geom.getMesh(), 11, 1);
DirichletBoundaryCondition(disp_field, "DISPLACEMENT_X", bc.begin(), bc.end(), 0.);
DirichletBoundaryCondition(disp_field, "DISPLACEMENT_Y", bc.begin(), bc.end(), 0.);
}
class MyNormal : public xEval<xVector> {
public:
void operator()(const xGeomElem* appro, const xGeomElem* integ, xVector& vec) const {
auto pt=integ->getXYZ();
double theta=atan2(pt(1), pt(0));
vec(0)=cos(theta);
vec(1)=sin(theta);
vec(2)=0.;
}
};
void QSFormulationSuperimposed::assembleNaturalEnv(LinearSystem& system) {
Formulation::assembleNaturalEnv(system);
NonUniformMaterialSensitivity_c<xTensor4> eval_stiffness("secant_matrix", variab_manager);
xFormLinearWithLoad<xGradOperator<xSymmetrize>, xEval<xTensor2> > form_linear(eval_stress_superimposed);
system.setAssemblerCoeff(-1.);
xEvalBinary<xMult<xTensor4, xTensor2, xTensor2>> eval_stress(eval_stiffness, eval_strain_superimposed);
xFormLinearWithLoad<xGradOperator<xSymmetrize>, xEval<xTensor2> > form_linear(eval_stress);
xFilteredRegion<xIter, xEntityFilter> superimposed(geom.begin(), geom.end(), xAcceptOnZoneID(101));
system.setAssemblerCoeff(-1.);
Assemble(form_linear, system.getAssembler(), geom.getIntegRuleSmart(), disp_field, superimposed.begin(), superimposed.end());
xEvalNormal eval_normal;
system.setAssemblerCoeff(1.);
MyNormal eval_normal;
xEvalBinary<xMult<xTensor2, xVector, xVector> > eval_traction(eval_stress_superimposed, eval_normal);
xFormLinearWithLoad<xValOperator<xIdentity<xVector> >, xEval<xVector> > lin(eval_traction);
xClassRegion bc(&geom.getMesh(), 12, 1);
Assemble(lin, system.getAssembler(), geom.getIntegRuleBasic(bc_integ_order), disp_field, bc.begin(), bc.end(), xUpperAdjacency());
system.setAssemblerCoeff(1.);
}
void QSFormulationSuperimposed::writeStrain(std::string name, std::string domain_name) {
......
......@@ -19,10 +19,6 @@ public:
const xfem::xEval<xfem::xTensor2>&, const xfem::xEval<xfem::xTensor2>&);
virtual ~QSFormulationSuperimposed();
// Constrains value by declaring it fixed (value state),
// gives the prescribed value from main.dat.
virtual void applyEssentialEnv(double);
// Assembles external force vector by reading main.dat file and
// assembles forces due to prescribed strain field.
virtual void assembleNaturalEnv(LinearSystem&);
......
......@@ -76,8 +76,6 @@ void ReadData(std::string main_file_name, xfem::xData& data, xParseData& parse_d
parse_data.registerString("disp_space_dim", "V2Dxy"); // selects the tensorial dimension of the space, see xSpaceFactory.h
parse_data.registerInt("disp_space_order", 1); // is the space order for the displacement field
parse_data.registerInt("disp_bc_integ_order", 1); // is the integ order used to integrate boundary conditions on displacement field
parse_data.registerInt("do_superposition", 0); // activate an additional contribution when computing strain in writeStrain
parse_data.registerTensor2("superimposed_strain", xTensor2(0.)); // additional strain
// mean field related
parse_data.registerString("mean_field_space_type", "PolyLagrange"); // selects phi and mean fields space (base space), see Util.h
parse_data.registerString("damage_shape", "Poly2Revert"); // selects damage shape function, see damageShapeFunctions.h
......@@ -87,7 +85,6 @@ void ReadData(std::string main_file_name, xfem::xData& data, xParseData& parse_d
// discretization related
parse_data.registerInt("do_crack_cut", 1); // cuts the crack and remeshes
parse_data.registerInt("do_bnd_in", 0); // use \Gamma_in instead of the interface between local and nonlocal (take domain bnd into account)
parse_data.registerInt("do_tls_solver_matrix_optim", 0); // the mass matrix involved in computation of mean fields is sometimes not needed
parse_data.registerInt("build_unconnected_parts", 0); // creates groups of unconnected domain parts, applies b.c. to all groups except first one
parse_data.registerDouble("fit_to_vertices_ratio", 1.e-2); // is a ratio to modulate the fit to vertices threshold value: epsilon=h*fit_to_vertices_ratio
parse_data.registerInt("remesh_aniso", 0); // anisotropic remesh
......@@ -96,7 +93,6 @@ void ReadData(std::string main_file_name, xfem::xData& data, xParseData& parse_d
parse_data.registerInt("nb_iter_max", 10); // is the max nb of iteration to compute equilibrium in a nonlinear context
parse_data.registerDouble("residual_norm_threshold", 1.e-4); // is the residual norm threshold in a nonlinear context
// export related
parse_data.registerString("curve_plot_file_ext", "none"); // selects the format that you want to plot your results
std::list<std::string> default_export_list={"damage", "1", "disp", "1"}; // for export_manager and export_sensor, see dedicated test case in Xtest
parse_data.registerListString("export_manager", default_export_list);
default_export_list.clear();
......
......@@ -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,53 @@ 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*dA/(A*A*A);
}
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*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 {
// d=min(d,0.99999);
return d/pow(sqrt(1.-d)+lda*pow(1.-sqrt(1.-d), 2.), 2.);