Commit 0d20941f authored by Kévin Moreau's avatar Kévin Moreau

This commit bring the following changes

* clean up unused feature.

it affects:
  - Export.cc
  - Import.cc

* improvement of transfer of displacement field after remeshing

it affects:
  - Formulation.cc

* correct a small bug that happens sometimes in restart, "mean" vectors
had wrong size after restart. It was always too long and filled with NaN.

it affects:
  - FormulationQS.cc

* had several "guards" that end computation and give an error message when
something goes wrong.

it affects:
  - TLSSolver.cc
  - FormulationQS.cc

* correct an error in the getGrad function of CohesiveLinearHardeningFunction

it affects:
  - MaterialFunction.cc

* swap orientation in comp mesh generation. It seems that it is test case dependent.
TO BE INVESTIGATED if your comp_mesh has two colours when plotting in GMSH.

it affects:
  - MeshGeneration.h

* correct a bug resulting from \Gamma_in recent commit. FastMarchingModeExtension
was updating phi based on \Gamma_in. The obtained phi was not anymore corresponding
to the iso-lc discretisation, it leads to damage equal to 1 inside the mechanical domain.
It leads to singular stiffness matrix.

it affects:
  - TLSSolver.cc
parent 90d823ac
......@@ -82,7 +82,6 @@ PostProcessing::PostProcessing(const xParseData& parse_data_, xEntityFilter filt
export_manager(GetNbDigits(parse_data.getInt("nb_step_max"))),
filter(filter) {
export_manager.init(parse_data, "export_manager");
createCurvesPlotFiles(parse_data.getString("curve_plot_file_ext"));
}
void PostProcessing::init(xMesh& mesh) {
......@@ -188,71 +187,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();
}
}
......@@ -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");
......
......@@ -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;
}
}
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)std::count_if(delta_phi.begin(), delta_phi.end(), [](double a){ return a>0.; }));
post_pro.exportOnTime("nb_active_modes", Observer::tell<int>("step"), (double)cnt);
}
void QSFormulation::updateLevelSetField() {
......@@ -426,11 +440,11 @@ 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;
// }
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;
......@@ -495,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());
}
......@@ -71,8 +71,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
......@@ -82,7 +80,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
......@@ -91,7 +88,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();
......
......@@ -227,7 +227,7 @@ double CohesiveLinearHardeningFunction::getVal(double d) const {
// 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;
return 1./(A*A)-2.*d*dA/(A*A*A);
}
double CohesiveLinearHardeningFunction::getGrad(double d) const {
......@@ -249,10 +249,11 @@ double CohesiveLinearHardeningFunction::getGrad(double d) const {
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 -4.*dA/(A*A*A)+6.*d*dA*dA/(A*A*A*A)-2.*d*d2A/(A*A*A);
}
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.);
}
......
......@@ -521,9 +521,9 @@ private:
// to propagate classification end
AOMD::mEntity* fd;
if(orientation(nds[0], nds[1], nds[2])) {
fd=comp_mesh.createFaceWithVertices((AOMD::mVertex*)nds[0], (AOMD::mVertex*)nds[1], (AOMD::mVertex*)nds[2], f_class);
} else {
fd=comp_mesh.createFaceWithVertices((AOMD::mVertex*)nds[1], (AOMD::mVertex*)nds[0], (AOMD::mVertex*)nds[2], f_class);
} else {
fd=comp_mesh.createFaceWithVertices((AOMD::mVertex*)nds[0], (AOMD::mVertex*)nds[1], (AOMD::mVertex*)nds[2], f_class);
}
fd->attachEntity(was_created_by_tag, fp);
sm_duplicated.add(fd);
......
......@@ -327,12 +327,28 @@ void TLSSolver::delocalizeMeanField() {
auto bnd_in_nonlocal=geom.getRegion("bnd_in_nonlocal");
auto int_in_nonlocal=geom.getRegion("int_in_nonlocal");
lalg::xDenseMatrix coeffs(bnd_in_nonlocal.size(0), nonlocal.size(0));
// TODO refactor this it is too complicated
std::function<void (mVertex* v, double val)> set_val;
if(epsilon_ratio==0) {
set_val=[](mVertex* v, double val){};
}
else {
set_val=[this](mVertex* v, double val){ this->field.setVal(v, val); };
}
FastMarchingModeExtension([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); },
set_val,
nonlocal,
bnd_in_nonlocal.begin(0), bnd_in_nonlocal.end(0),
int_in_nonlocal.begin(0), int_in_nonlocal.end(0),
coeffs, epsilon_ratio);
const auto* vals=coeffs.ReturnValPointer();
const double min=*std::min_element(vals, vals+coeffs.nline()*coeffs.ncolumn());
if(min<0.) {
std::cout<<"Error: coefficients in nonlocality_matrix must be positive!"<<std::endl;
std::cout<<" There is maybe a bug in fast marching or in the input data"<<std::endl;
std::abort();
}
geom.deleteBndIntInNonlocal();
if(nonlocality_matrix) {
......
......@@ -36,7 +36,8 @@ xfem::xSpaceFactoryBase* SpaceTypeSelector(std::string name) {
return new xfem::xSpaceFactory<xfem::xSpacePolynomialLagrange, STT>();
}
else {
std::cout<<"Error, no space type named "<<name<<std::endl;
std::cout<<"Error: no space type named "<<name<<std::endl;
std::abort();
}
}
......
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