Commit 5bed31d2 authored by Kevin Moreau's avatar Kevin Moreau

Correct a bug in restart in FormulationQS.cc.

A line of code was missing for a restart after the emergence of a crack.
Material variables are now updated after element duplication.
Refactorise and comment a little more.



git-svn-id: https://svn.ec-nantes.fr/eXlibris/Applis/TLSDuctile@2347 fbbead7c-fb4d-4173-aa67-51132c73c120
parent ca797a50
......@@ -20,7 +20,6 @@
#include "xCSRVector.h"
#include "xCSRMatrix.h"
// #include "xDenseMatrix.h"
#include "MaterialCommand.h"
#include "NonUniformMaterialSensitivity.h"
......@@ -71,22 +70,7 @@ QSFormulation::QSFormulation(TLSGeom& geom, TLSSolver& tls_solver,
declareDofs();
if(pre_pro.doRestart()) {
pre_pro.loadField("disp", geom.getMesh(), disp_field);
xEvalField<xIdentity<xVector> > eval_disp(disp_field);
post_pro.exportOnSpace("restored_disp_basic", pre_pro.getStep(), eval_disp, geom.getIntegRuleBasic(0), geom.begin(), geom.end());
tls_solver.swapOldAndCurrentLevelSetField();
updateMaterialVariables(tls_solver.getDamageEval());
writeStrain("eps_ref");
tls_solver.swapOldAndCurrentLevelSetField();
updateMaterialVariables(tls_solver.getDamageEval());
if(tls_solver.isCrackable()) {
geom.buildFullyDamagedSupport(tls_solver.getPhiMinusLcEval());
tls_solver.cutCrack();
tls_solver.transferLevelSetField();
transferDispField();
transferMaterialVariables();
}
Observer::listen("step", pre_pro.getStep());
restoreStep(pre_pro);
}
if(is_nonlinear) {
......@@ -122,6 +106,9 @@ void QSFormulation::updateCurrentState() {
computeDeltaLevelSetField();
updateLevelSetField();
updateMaterialVariables(tls_solver.getDamageEval());
if(is_nonlinear) {
geom.buildLinearAndNonLinear(tls_solver.getPhiEval());
}
is_forced=false;
Display::stepOut();
}
......@@ -139,6 +126,9 @@ bool QSFormulation::updateCurrentStateForced() {
tls_solver.writeOldToCurrentLevelSetField();
updateLevelSetFieldForced();
updateMaterialVariables(tls_solver.getDamageEval());
if(is_nonlinear) {
geom.buildLinearAndNonLinear(tls_solver.getPhiEval());
}
is_forced=true;
Display::stepOut();
return true;
......@@ -221,6 +211,8 @@ void QSFormulation::computeLSAdvanceMeanFields() {
xEvalBinary<std::multiplies<double> > eval_beta_dd_over_dphi(eval_beta, eval_dd_over_dphi);
tls_solver.updateMeanField("mean_alpha", eval_dd_over_dphi, eval_alpha_dd_over_dphi);
tls_solver.updateMeanField("mean_beta", eval_dd_over_dphi, eval_beta_dd_over_dphi);
// tls_solver.exportMeanField("mean_alpha", Observer::tell<int>("step"));
// tls_solver.exportMeanField("mean_beta", Observer::tell<int>("step"));
}
void QSFormulation::computeEquilibriumMeanFields() {
......@@ -230,6 +222,7 @@ void QSFormulation::computeEquilibriumMeanFields() {
NonUniformMaterialSensitivity_c<double> eval_str(str, variab_manager);
xEvalBinary<std::multiplies<double> > eval_str_dd_over_dphi(eval_str, eval_dd_over_dphi);
tls_solver.updateMeanField("mean_"+str, eval_dd_over_dphi, eval_str_dd_over_dphi);
// tls_solver.exportMeanField("mean_"+str, Observer::tell<int>("step"));
}
}
......@@ -253,45 +246,32 @@ void QSFormulation::updateLevelSetField() {
tls_solver.updateLevelSetField();
tls_solver.delocalizeLevelSetField();
if(tls_solver.isDelocalized()) {
tls_solver.reinitLevelSetField();
if(tls_solver.isCrackable()) {
clearDispMeasField();
geom.buildFullyDamagedSupport(tls_solver.getPhiMinusLcEval());
deleteMaterialVariables("fds");
tls_solver.cutCrack();
tls_solver.transferLevelSetField();
transferDispField();
transferMaterialVariables();
declareDispMeasField();
}
tls_solver.delocalizeMeanField();
}
if(is_nonlinear) {
geom.buildLinearAndNonLinear(tls_solver.getPhiEval());
reinitLevelSetField();
}
}
void QSFormulation::updateLevelSetFieldForced() {
tls_solver.delocalizeLevelSetFieldForced(delta_phi_max);
if(tls_solver.isNewlyDelocalized()) {
tls_solver.reinitLevelSetField();
if(tls_solver.isCrackable()) {
clearDispMeasField();
geom.buildFullyDamagedSupport(tls_solver.getPhiMinusLcEval());
deleteMaterialVariables("fds");
tls_solver.cutCrack();
tls_solver.transferLevelSetField();
transferDispField();
transferMaterialVariables();
declareDispMeasField();
}
tls_solver.delocalizeMeanField();
}
if(is_nonlinear) {
geom.buildLinearAndNonLinear(tls_solver.getPhiEval());
reinitLevelSetField();
}
}
void QSFormulation::reinitLevelSetField() {
tls_solver.reinitLevelSetField();
if(tls_solver.isCrackable()) {
clearDispMeasField();
geom.buildFullyDamagedSupport(tls_solver.getPhiMinusLcEval());
deleteMaterialVariables("fds");
tls_solver.cutCrack();
tls_solver.transferLevelSetField();
transferDispField();
transferMaterialVariables();
declareDispMeasField();
}
tls_solver.delocalizeMeanField();
}
void QSFormulation::writeDisp(const xCSRVector& vec, std::string name) {
if(name=="disp_ref") {
CustomValueDoubleSelector::which_val=0;
......@@ -440,36 +420,6 @@ bool QSFormulation::exportStep(int step) {
tls_solver.exportStep(step);
// error
/*
const xCSRVector& a=tls_solver.getMeanVector("mean_a");
const xCSRVector& b=tls_solver.getMeanVector("mean_b");
const xCSRVector& c=tls_solver.getMeanVector("mean_c");
const xCSRVector& delta_phi=tls_solver.getMeanVector("delta_phi");
const int nb_mode=delta_phi.size();
xDenseMatrix mass(nb_mode, nb_mode);
xAssemblerBasic<lalg::xDenseMatrix> assembler(mass);
tls_solver.assembleMassMatrix(assembler);
xCSRVector f(nb_mode);
for(int i=0; i<nb_mode; ++i) {
f[i]=std::fabs(a[i]*load_factor*load_factor+b[i]*load_factor+c[i]);
}
xCSRVector tmp(nb_mode);
const int INCX = 1;
const int INCY = 1;
const double alpha = 1.;
const double beta = 0.;
dgemv_("N", &nb_mode, &nb_mode, &alpha, const_cast<double*>(mass.ReturnValPointer()),
&nb_mode, &delta_phi[0], &INCX, &beta, &tmp[0], &INCY);
xCSRVector tmp2(nb_mode);
xCSRVector one(nb_mode);
std::fill(one.begin(), one.end(), 1.);
dgemv_("N", &nb_mode, &nb_mode, &alpha, const_cast<double*>(mass.ReturnValPointer()),
&nb_mode, &one[0], &INCX, &beta, &tmp2[0], &INCY);
const double error=dot(f, tmp)/dot(one, tmp2);
post_pro.exportOnTime("error", (double)step, error);
*/
if(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;
......@@ -515,3 +465,34 @@ void QSFormulation::saveStep(int step) {
tls_solver.saveStep(step);
post_pro.archiveSave(step);
}
void QSFormulation::restoreStep(PreProcessing& pre_pro) {
Observer::listen("step", pre_pro.getStep());
pre_pro.loadField("disp", geom.getMesh(), disp_field);
xEvalField<xIdentity<xVector> > eval_disp(disp_field);
post_pro.exportOnSpace("restored_disp_basic", pre_pro.getStep(), eval_disp, geom.getIntegRuleBasic(0), geom.begin(), geom.end());
tls_solver.swapOldAndCurrentLevelSetField();
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");
}
tls_solver.delocalizeMeanField();
}
......@@ -17,35 +17,40 @@ public:
PreProcessing&, PostProcessing&);
virtual ~QSFormulation();
// Updates level set field and material variables
// Updates level set field, and material variables. In formulations,
// update is to take in a large sense (it also delocalizes,
// reinitializes, enrichs, transfers data).
virtual void updateCurrentState();
// Perturbs level set field and updates material variables
// Perturbs level set field and updates material variables.
virtual bool updateCurrentStateForced();
// Computes load factor variation
// Computes load factor variation.
virtual void computeLoadFactorVariation();
// Computes load factor
// Computes load factor.
double computeLoadFactor();
// Writes to disp field dof values from xCSRVector, a specific
// disp field is selected by the string
// disp field is selected by the string.
void writeDisp(const lalg::xCSRVector&, std::string);
// Appends to disp field dof values from xCSRVector.
void addDisp(const lalg::xCSRVector&, std::string);
// Writes current values of level set field and material variables
// to old values
// to old values.
void writeCurrentToOld();
// Exports for post processing
// Exports for post processing.
virtual bool exportStep(int);
// Saves for restoration
// Saves for restoration.
virtual void saveStep(int);
protected:
// Restores phi field and old phi field, computes delocalization matrix.
void restoreStep(PreProcessing&);
// Calls tls solver to compute all required mean fields to update phi.
// Required mean fields are declared in the constructor.
void computeLSAdvanceMeanFields();
......@@ -57,13 +62,25 @@ protected:
// Computes delta phi and writes it to tls solver.
virtual void computeDeltaLevelSetField();
void updateLevelSetFieldForced();
// Updates, delocalizes and reinitializes phi.
virtual void updateLevelSetField();
// Updates at one node only (Perturbs), delocalizes and reinitializes phi.
void updateLevelSetFieldForced();
// Reinitializes phi.
void reinitLevelSetField();
// Clears measurement displacement field.
void clearDispMeasField();
// Declares measurement displacement field.
void declareDispMeasField();
// Computes dual (energetically) force of the imposed displacement.
xfem::xVector computeRefForce();
// Computes dual (energetically) displacement of the imposed stress.
xfem::xVector computeRefDisp();
protected:
TLSSolver& tls_solver;
......
......@@ -152,10 +152,11 @@ TLSGeom::TLSGeom(xData& data,
xMesh* comp_mesh;
std::list<mEntity*> nonlocal_list;
if(pre_pro.doRestart()) {
// to do a clean restart without loss uncomment this 2 lines (only if there is no crack)
// data.ReadMesh();
// comp_mesh=data.mesh;// TODO stored mesh as node position perturbed w.r.t. original mesh
// comp_mesh=data.mesh;
comp_mesh=new xMesh;
pre_pro.loadMesh("comp", *comp_mesh);
pre_pro.loadMesh("comp", *comp_mesh);// TODO, Warning: the loaded mesh has node location slightly different from original one. This is because we export the mesh in saveStep in ASCII.
pre_pro.loadDomain("nonlocal", *comp_mesh, nonlocal_list);
}
else {
......@@ -800,11 +801,18 @@ void TLSGeom::buildCompMesh() {
comp_mesh.modifyAllState();// TODO maybe improvable
comp_mesh.bdryLinkSetup();// TODO maybe improvable
for(auto it=sm_duplicated.begin(0); it!=sm_duplicated.end(0); ++it) {
AddNeighbors(*it, sm_duplicated, 1); // needed for fast marching
}
xRegion integ(&getMesh("integ"));
std::for_each(integ.begin(), integ.end(), PartitionCreator(partition_tag, xAttached(duplicated_tag)));
auto step=Observer::tell<int>("step");
post_pro.exportOnSpace("comp_mesh", step, comp_mesh);
xEvalConstant<double> eval_one(1.);
xIntegrationRuleBasic integ_rule(0);
post_pro.exportOnSpace("sm_duplicated", step, eval_one, integ_rule, sm_duplicated.begin(), sm_duplicated.end());
}
void TLSGeom::createIntegrationPoints() {
......@@ -851,96 +859,93 @@ void TLSGeom::cleanCrack() {
}
void TLSGeom::saveStep(int step) {
if(isRegisteredMesh("integ")) {
post_pro.saveMesh("integ", step, getMesh("integ"));
}
xMesh& mesh=getMesh();
post_pro.saveMesh("comp", step, mesh);
xRegion nonlocal(&mesh.getSubMesh("nonlocal"));
post_pro.saveDomain("nonlocal", step, nonlocal.begin(), nonlocal.end());
}
void classifyUnconnectedRegions(const xMesh& m, const xSubMesh& sm, int dim, std::list<xSubMesh*>& sm_disc_list) {
assert(sm_disc_list.empty());
xAcceptInSubMesh is_in_sm(sm);
if(dim==0) {
dim=m.dim();
for(auto it=sm.begin(0); it!=sm.end(0); ++it) {
mEntity* n=*it;
// look for new region
bool is_new_region=true;
for(std::list<xSubMesh*>::iterator itl=sm_disc_list.begin(); itl!=sm_disc_list.end(); ++itl) {
if((*itl)->find(n)) {
is_new_region=false;
break;
}
}
// add new region and fill it
if(is_new_region) {
std::ostringstream oss;
oss<<sm_disc_list.size();
sm_disc_list.push_back(new xSubMesh(sm.getName()+"_disc_"+oss.str(), m));
xSubMesh* sm_disc=sm_disc_list.back();
xAcceptInSubMesh is_in_disc(*sm_disc);
std::list<mEntity*> totreat;
totreat.push_back(n);
sm_disc->add(n);
while(!totreat.empty()) {
mEntity* nn=totreat.back();
totreat.pop_back();
for(int i=0; i<nn->size(dim); ++i) {
mEntity* f=nn->get(dim,i);
if(is_in_sm(f)) {
for(int j=0; j<f->size(0); ++j) {
mEntity* nnn=f->get(0,j);
if(nnn!=nn && is_in_sm(nnn) && !is_in_disc(nnn)) {
totreat.push_back(nnn);
sm_disc->add(nnn);
}
}
}
}
}
}
}
}
else {
assert(dim==2||dim==3);
for(auto it=sm.begin(dim); it!=sm.end(dim); ++it) {
mEntity* f=*it;
// look for new region
bool is_new_region=true;
for(std::list<xSubMesh*>::iterator itl=sm_disc_list.begin(); itl!=sm_disc_list.end(); ++itl) {
if((*itl)->find(f)) {
is_new_region=false;
break;
}
}
// add new region and fill it
if(is_new_region) {
std::ostringstream oss;
oss<<sm_disc_list.size();
sm_disc_list.push_back(new xSubMesh(sm.getName()+"_disc_"+oss.str(), m));
xSubMesh* sm_disc=sm_disc_list.back();
xAcceptInSubMesh is_in_disc(*sm_disc);
std::list<mEntity*> totreat;
totreat.push_back(f);
sm_disc->add(f);
while(!totreat.empty()) {
mEntity* ff=totreat.back();
totreat.pop_back();
for(int i=0; i<ff->size(0); ++i) {
mEntity* n=ff->get(0,i);
for(int j=0; j<n->size(dim); ++j) {
mEntity* fff=n->get(dim,j);
if(fff!=ff && is_in_sm(fff) && !is_in_disc(fff)) {
totreat.push_back(fff);
sm_disc->add(fff);
}
}
}
}
}
}
}
}
// void classifyUnconnectedRegions(const xMesh& m, const xSubMesh& sm, int dim, std::list<xSubMesh*>& sm_disc_list) {
// assert(sm_disc_list.empty());
// xAcceptInSubMesh is_in_sm(sm);
// if(dim==0) {
// dim=m.dim();
// for(auto it=sm.begin(0); it!=sm.end(0); ++it) {
// mEntity* n=*it;
// // look for new region
// bool is_new_region=true;
// for(std::list<xSubMesh*>::iterator itl=sm_disc_list.begin(); itl!=sm_disc_list.end(); ++itl) {
// if((*itl)->find(n)) {
// is_new_region=false;
// break;
// }
// }
// // add new region and fill it
// if(is_new_region) {
// std::ostringstream oss;
// oss<<sm_disc_list.size();
// sm_disc_list.push_back(new xSubMesh(sm.getName()+"_disc_"+oss.str(), m));
// xSubMesh* sm_disc=sm_disc_list.back();
// xAcceptInSubMesh is_in_disc(*sm_disc);
// std::list<mEntity*> totreat;
// totreat.push_back(n);
// sm_disc->add(n);
// while(!totreat.empty()) {
// mEntity* nn=totreat.back();
// totreat.pop_back();
// for(int i=0; i<nn->size(dim); ++i) {
// mEntity* f=nn->get(dim,i);
// if(is_in_sm(f)) {
// for(int j=0; j<f->size(0); ++j) {
// mEntity* nnn=f->get(0,j);
// if(nnn!=nn && is_in_sm(nnn) && !is_in_disc(nnn)) {
// totreat.push_back(nnn);
// sm_disc->add(nnn);
// }
// }
// }
// }
// }
// }
// }
// }
// else {
// assert(dim==2||dim==3);
// for(auto it=sm.begin(dim); it!=sm.end(dim); ++it) {
// mEntity* f=*it;
// // look for new region
// bool is_new_region=true;
// for(std::list<xSubMesh*>::iterator itl=sm_disc_list.begin(); itl!=sm_disc_list.end(); ++itl) {
// if((*itl)->find(f)) {
// is_new_region=false;
// break;
// }
// }
// // add new region and fill it
// if(is_new_region) {
// std::ostringstream oss;
// oss<<sm_disc_list.size();
// sm_disc_list.push_back(new xSubMesh(sm.getName()+"_disc_"+oss.str(), m));
// xSubMesh* sm_disc=sm_disc_list.back();
// xAcceptInSubMesh is_in_disc(*sm_disc);
// std::list<mEntity*> totreat;
// totreat.push_back(f);
// sm_disc->add(f);
// while(!totreat.empty()) {
// mEntity* ff=totreat.back();
// totreat.pop_back();
// for(int i=0; i<ff->size(0); ++i) {
// mEntity* n=ff->get(0,i);
// for(int j=0; j<n->size(dim); ++j) {
// mEntity* fff=n->get(dim,j);
// if(fff!=ff && is_in_sm(fff) && !is_in_disc(fff)) {
// totreat.push_back(fff);
// sm_disc->add(fff);
// }
// }
// }
// }
// }
// }
// }
// }
......@@ -80,26 +80,6 @@ TLSSolver::~TLSSolver() {
delete space_factory;
}
void TLSSolver::restoreStep(PreProcessing& pre_pro) {
Observer::listen("step", pre_pro.getStep());
auto integ_rule=geom.getIntegRuleBasic(0);
auto begin=geom.begin("tls");
auto end=geom.end("tls");
OldAndCurrent_c::old();
pre_pro.loadField("old_phi", geom.getMesh(), field);
post_pro.exportOnSpace("restored_old_phi", pre_pro.getStep(), eval_phi, integ_rule, begin, end);
OldAndCurrent_c::current();
pre_pro.loadField("phi", geom.getMesh(), field);
post_pro.exportOnSpace("restored_phi", pre_pro.getStep(), eval_phi, integ_rule, begin, end);
xEvalGradField<xIdentity<xVector> > eval_grad_phi(field);
geom.updateBndAndInt(eval_grad_phi);
if(isDelocalized()) {
computeModesFastMarching();
}
}
void TLSSolver::deleteState() {
DeleteState(field, geom.begin(0, "tls"), geom.end(0, "tls"));
double_manager.clear_subset(dofs);
......@@ -296,9 +276,15 @@ void TLSSolver::transferLevelSetField() {
transferDomains();
declareField("duplicated");
deleteState();
TransferField(field, field,
[](xValue<double>* v1, xValue<double>* v2) { v2->setVal(v1->getVal()); },
geom.begin(0, "duplicated"), geom.end(0, "duplicated"));
// TransferField(field, field,
// [](xValue<double>* v1, xValue<double>* v2) { v2->setVal(v1->getVal()); },
// geom.begin(0, "duplicated"), geom.end(0, "duplicated"));
// it is better than a transfer since new duplicated nodes (same location) receive different values (one from left, one from right in 1D)
FastMarchingReinit(field,
geom.getRegion("duplicated"),
geom.begin(0, "bnd_fds"), geom.end(0, "bnd_fds"),
geom.begin(0, "duplicated"), geom.end(0, "duplicated"),
[](const mEntity& v) { return 1.; }, false);
deleteField("int_fds");
declareState();
post_pro.exportOnSpace("transfered_phi", Observer::tell<int>("step"), eval_phi, geom.getIntegRuleBasic(0), geom.begin("tls"), geom.end("tls"));
......@@ -396,13 +382,6 @@ bool TLSSolver::isDelocalized() const {
return isNewlyDelocalized() || geom.isDelocalized();
}
void TLSSolver::saveStep(int step) {
OldAndCurrent_c::old();
post_pro.saveField("old_phi", step, field, geom.getRegion("tls"));
OldAndCurrent_c::current();
post_pro.saveField("phi", step, field, geom.getRegion("tls"));
}
void TLSSolver::exportStep(int step) {
auto integ_rule=geom.getIntegRulePartition(0);
auto begin=geom.begin("tls");
......@@ -437,7 +416,26 @@ const xEval<double>& TLSSolver::getPhiMinusLcEval() const {
return eval_phi_minus_lc;
}
// void TLSSolver::assembleMassMatrix(xAssemblerBasic<lalg::xDenseMatrix>& assembler) {
// xFormBilinearWithoutLaw<xValOperator<xIdentity<double>>, xValOperator<xIdentity<double>>> form_bilinear;
// Assemble(form_bilinear, assembler, geom.getIntegRuleBasic(), mean_field, mean_field, geom.begin("tls"), geom.end("tls"));
// }
void TLSSolver::saveStep(int step) {
OldAndCurrent_c::old();
post_pro.saveField("old_phi", step, field, geom.getRegion("tls"));
OldAndCurrent_c::current();
post_pro.saveField("phi", step, field, geom.getRegion("tls"));
}
void TLSSolver::restoreStep(PreProcessing& pre_pro) {
Observer::listen("step", pre_pro.getStep());
auto integ_rule=geom.getIntegRuleBasic(0);
auto begin=geom.begin("tls");
auto end=geom.end("tls");
OldAndCurrent_c::old();
pre_pro.loadField("old_phi", geom.getMesh(), field);
post_pro.exportOnSpace("restored_old_phi", pre_pro.getStep(), eval_phi, integ_rule, begin, end);
OldAndCurrent_c::current();
pre_pro.loadField("phi", geom.getMesh(), field);
post_pro.exportOnSpace("restored_phi", pre_pro.getStep(), eval_phi, integ_rule, begin, end);
xEvalGradField<xIdentity<xVector> > eval_grad_phi(field);
geom.updateBndAndInt(eval_grad_phi);
}
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