Commit 4808c191 authored by Kevin Moreau's avatar Kevin Moreau

Mostly refactoring.

Changed the transfer of level set field by a reinit level set field
on duplicated nodes.



git-svn-id: https://svn.ec-nantes.fr/eXlibris/Applis/TLSDuctile@2361 fbbead7c-fb4d-4173-aa67-51132c73c120
parent 4d372035
......@@ -12,13 +12,46 @@
#include "xDenseMatrix.h"
//
#ifdef HAVE_FASTMARCHING
#include "gmshinputoutput.h"
// #include "gmshinputoutput.h"
#include "meshinterfacexRegion.h"
#include "FastMarching.h"
// #include "FM.h"
#endif
// std
#include <functional>
class modesValues{
public:
modesValues(){};
modesValues(size_t n, double val):values(n,val){};
modesValues(const modesValues &in):values(in.values){};
double & operator[](const size_t &i){ return values[i];}
modesValues & operator*=(const double &scal){
for_each(values.begin(), values.end(), [&scal](double & val){ val*= scal;} );
return *this;
}
modesValues & operator+=(const modesValues &in){
transform( values.begin(), values.end(), in.values.begin(), values.begin(), [] (double val1, double val2 ){return val1+val2;});
return *this;
}
modesValues& operator=(const modesValues &in) {
values=in.values;
return *this;
}
friend modesValues operator+(modesValues lhs, const modesValues& rhs)
{
return lhs += rhs; // reuse compound assignment and return the result by value
}
friend modesValues operator*(const double &lhs, modesValues rhs){
return rhs *= lhs;
}
friend modesValues operator*(modesValues rhs, const double &lhs){
return rhs *= lhs;
}
std::vector<double> values;
private :
};
template<typename FIELD, typename ITERNODE>
void FastMarchingReinit(FIELD& phi,
const xfem::xRegion& region,
......@@ -50,17 +83,23 @@ void FastMarchingReinit(FIELD& phi,
}
meshinterfacexRegion mi(region);
entitystorage<meshinterfacexRegion, AOMD::mVertex, double> ls(mi);
// entitystorage<meshinterfacexRegion, AOMD::mVertex, vector3d<double> > gls(mi);
// entitystorage<meshinterfacexRegion, AOMD::mVertex, int> vn(mi);
for(ITERNODE it=known_begin; it!=known_end; ++it) {
auto v=*it;
std::vector<double> vals;
phi.getVals(v, vals);
ls.set(*static_cast<AOMD::mVertex*>(v), sign*vals[0]+max-min);
}
// for(ITERNODE it=tofind_begin; it!=tofind_end; ++it) {
// ls.set(*static_cast<AOMD::mVertex*>(*it), 1000.);
// }
entitytovertexiteratorconvertor<ITERNODE> known_conv_begin(known_begin);
entitytovertexiteratorconvertor<ITERNODE> known_conv_end(known_end);
entitytovertexiteratorconvertor<ITERNODE> tofind_conv_begin(tofind_begin);
entitytovertexiteratorconvertor<ITERNODE> tofind_conv_end(tofind_end);
fmeik(mi, ls, known_conv_begin, known_conv_end, tofind_conv_begin, tofind_conv_end, f_func);
// fmeik(mi, ls, known_conv_begin, known_conv_end, f_func, gls, vn);
for(ITERNODE it=tofind_begin; it!=tofind_end; ++it) {
auto v=*it;
double val;
......@@ -70,38 +109,6 @@ void FastMarchingReinit(FIELD& phi,
#endif
}
class modesValues{
public:
modesValues(){};
modesValues(size_t n, double val):values(n,val){};
modesValues(const modesValues &in):values(in.values){};
double & operator[](const size_t &i){ return values[i];}
modesValues & operator*=(const double &scal){
for_each(values.begin(), values.end(), [&scal](double & val){ val*= scal;} );
return *this;
}
modesValues & operator+=(const modesValues &in){
transform( values.begin(), values.end(), in.values.begin(), values.begin(), [] (double val1, double val2 ){return val1+val2;});
return *this;
}
modesValues& operator=(const modesValues &in) {
values=in.values;
return *this;
}
friend modesValues operator+(modesValues lhs, const modesValues& rhs)
{
return lhs += rhs; // reuse compound assignment and return the result by value
}
friend modesValues operator*(const double &lhs, modesValues rhs){
return rhs *= lhs;
}
friend modesValues operator*(modesValues rhs, const double &lhs){
return rhs *= lhs;
}
std::vector<double> values;
private :
};
template<typename FIELD, typename ITERNODE>
void FastMarchingModeExtension(FIELD& phi,
const xfem::xRegion& region,
......@@ -114,6 +121,7 @@ void FastMarchingModeExtension(FIELD& phi,
#else
meshinterfacexRegion mi(region);
entitystorage<meshinterfacexRegion, AOMD::mVertex, double> ls(mi);
// entitystorage<meshinterfacexRegion, AOMD::mVertex, vector3d<double> > gls(mi);
entitystorage<meshinterfacexRegion, AOMD::mVertex, modesValues> vn(mi);
const int nb_modes=std::distance(known_begin, known_end);
int i=0;
......@@ -126,12 +134,18 @@ void FastMarchingModeExtension(FIELD& phi,
mv[i]=1.;
vn.set(*static_cast<AOMD::mVertex*>(v), mv);
}
// for(ITERNODE it=tofind_begin; it!=tofind_end; ++it, ++i) {
// ls.set(*static_cast<AOMD::mVertex*>(*it), 1000.);
// modesValues mv(nb_modes, 0.);
// vn.set(*static_cast<AOMD::mVertex*>(*it), mv);
// }
entitytovertexiteratorconvertor<ITERNODE> known_conv_begin(known_begin);
entitytovertexiteratorconvertor<ITERNODE> known_conv_end(known_end);
entitytovertexiteratorconvertor<ITERNODE> tofind_conv_begin(tofind_begin);
entitytovertexiteratorconvertor<ITERNODE> tofind_conv_end(tofind_end);
std::function<double (const AOMD::mVertex&)> f_func=[](const AOMD::mVertex& v){return 1.;};
fmeik(mi, ls, known_conv_begin, known_conv_end, tofind_conv_begin, tofind_conv_end, f_func, vn);
// fmeik(mi, ls, known_conv_begin, known_conv_end, f_func, gls, vn);
coeffs.resize(nb_modes, region.size(0));
int j=0;
for(auto it=region.begin(0); it!=region.end(0); ++it, ++j) {
......
......@@ -70,7 +70,7 @@ QSFormulation::QSFormulation(TLSGeom& geom, TLSSolver& tls_solver,
std::cout<<"delta phi max "<<delta_phi_max<<std::endl;
}
QSFormulation::~QSFormulation(){}
QSFormulation::~QSFormulation() {}
void QSFormulation::initializeValueDouble() {
std::string mat_class=getMaterialClass();
......@@ -279,6 +279,7 @@ void QSFormulation::reinitLevelSetField() {
tls_solver.transferLevelSetField();
transferDispField();
transferMaterialVariables();
geom.deleteFullyDamagedSupport();
declareDispMeasField();
}
tls_solver.delocalizeMeanField();
......@@ -477,8 +478,6 @@ void QSFormulation::saveStep(int 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());
......
......@@ -130,6 +130,7 @@ PreProcessing::PreProcessing(std::string archive_filename) :
iss.read((char*)&old_load_factor, sizeof(double));
iss.read((char*)&load_factor, sizeof(double));
iss.read((char*)&load_factor_variation, sizeof(double));
Observer::listen("step", restart_step);
Observer::listen("old_load_factor", old_load_factor);
Observer::listen("load_factor", load_factor);
Observer::listen("load_factor_variation", load_factor_variation);
......
......@@ -320,17 +320,17 @@ bool LessThanId::operator()(int id, mEntity* fi) const {
NodeDispatcher::NodeDispatcher(const xMesh& comp_mesh,
const unsigned int was_created_by_tag,
const unsigned int is_crack_tag,
const unsigned int partition_tag) :
const unsigned int partition_tag,
const unsigned int duplicated_node_tag) :
dim(comp_mesh.dim()),
was_created_by_tag(was_created_by_tag),
is_crack_tag(is_crack_tag),
partition_tag(partition_tag),
support_tag(AOMD_Util::Instance()->lookupMeshDataId("support")),
duplicated_node_tag(AOMD_Util::Instance()->lookupMeshDataId("duplicated_node"))
duplicated_node_tag(duplicated_node_tag),
support_tag(AOMD_Util::Instance()->newMeshDataId("support"))
{}
NodeDispatcher::~NodeDispatcher() {
// AOMD_Util::Instance()->deleteMeshDataId(duplicated_node_tag);
AOMD_Util::Instance()->deleteMeshDataId(support_tag);
}
......@@ -368,20 +368,19 @@ 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 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),
partition_tag(partition_tag),
support_tag(AOMD_Util::Instance()->lookupMeshDataId("support")),
duplicated_node_tag(AOMD_Util::Instance()->lookupMeshDataId("duplicated_node"))
duplicated_node_tag(duplicated_node_tag),
support_tag(AOMD_Util::Instance()->newMeshDataId("support"))
{}
NodeDuplicator::~NodeDuplicator() {
// AOMD_Util::Instance()->deleteMeshDataId(duplicated_node_tag);
AOMD_Util::Instance()->deleteMeshDataId(support_tag);
}
......@@ -420,19 +419,19 @@ FaceDuplicator::FaceDuplicator(xMesh& comp_mesh,
const unsigned int was_created_by_tag,
const unsigned int is_crack_tag,
const unsigned int partition_tag,
const unsigned int duplicated_tag) :
const unsigned int duplicated_tag,
const unsigned int duplicated_node_tag) :
comp_mesh(comp_mesh), sm_duplicated(sm_duplicated), id(0),
dim(comp_mesh.dim()),
was_created_by_tag(was_created_by_tag),
is_crack_tag(is_crack_tag),
partition_tag(partition_tag),
duplicated_tag(duplicated_tag),
local_id_tag(AOMD_Util::Instance()->lookupMeshDataId("local_id")),
duplicated_node_tag(AOMD_Util::Instance()->lookupMeshDataId("duplicated_node")) {} // attention detruire et recreer est dangeureux
duplicated_node_tag(duplicated_node_tag),
local_id_tag(AOMD_Util::Instance()->newMeshDataId("local_id")) {}
FaceDuplicator::~FaceDuplicator() {
AOMD_Util::Instance()->deleteMeshDataId(local_id_tag);
AOMD_Util::Instance()->deleteMeshDataId(duplicated_node_tag);
}
void FaceDuplicator::operator()(mEntity* fp) {
......
......@@ -385,6 +385,7 @@ struct LessThanId {
class NodeDispatcher {
public:
NodeDispatcher(const xfem::xMesh&,
const unsigned int,
const unsigned int,
const unsigned int,
const unsigned int);
......@@ -411,6 +412,7 @@ public:
xfem::xSubMesh&,
const unsigned int,
const unsigned int,
const unsigned int,
const unsigned int);
~NodeDuplicator();
void operator()(AOMD::mEntity*);
......@@ -441,6 +443,7 @@ public:
const unsigned int,
const unsigned int,
const unsigned int,
const unsigned int,
const unsigned int);
~FaceDuplicator();
void operator()(AOMD::mEntity*);
......
This diff is collapsed.
......@@ -61,7 +61,9 @@ public:
void updateDomains(std::list<AOMD::mEntity*>&);
void transferDomains();
void updateBndAndInt(const xfem::xEval<xfem::xVector>&);
void deleteBndAndInt();
void buildFullyDamagedSupport(const xfem::xEval<double>&);
void deleteFullyDamagedSupport();
void buildLinearAndNonLinear(const xfem::xEval<double>&);
bool isDelocalized() const;
......@@ -121,9 +123,9 @@ private:
double h, body_characteristic_length;
const unsigned int was_created_by_tag;
const unsigned int duplicated_tag;
const unsigned int partition_tag;
const unsigned int duplicated_partition_tag;
const unsigned int duplicated_tag;
const unsigned int duplicated_node_tag;
const unsigned int is_crack_tag;
};
......
......@@ -54,13 +54,16 @@ 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"))) {
space_factory->setSpaceProductOrder(1);
field.insert(space_factory->getSpace("fem"));
declareField();
declareState();
if(pre_pro.doRestart()) {
restoreStep(pre_pro);
}
registerMeanVector("delta_phi");
#ifndef TLSMATRIXOPTIM
registerMeanVector("diag");
......@@ -79,26 +82,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);
......@@ -174,13 +157,6 @@ void TLSSolver::swapOldAndCurrentLevelSetField() {
void TLSSolver::updateLevelSetField() {
ConvertMeanVectorModalToNodal(nonlocality_matrix, getMeanVector("delta_phi"), buf_nodal_vec);
Visit(xAddSolutionVisitor<>(buf_nodal_vec.begin()), double_manager.begin(dofs), double_manager.end(dofs));
auto integ_rule=geom.getIntegRuleBasic(0);
auto begin=geom.begin("tls");
auto end=geom.end("tls");
auto step=Observer::tell<int>("step");
post_pro.exportOnSpace("phi_after_update", step, eval_phi, integ_rule, begin, end);
xEvalGradField<xNorm> eval_norm_grad_phi(field);
post_pro.exportOnSpace("norm_grad_phi_after_update", step, eval_norm_grad_phi, integ_rule, begin, end);
#ifndef TLSMATRIXOPTIM
assemble_mat=true;
#endif
......@@ -204,9 +180,6 @@ void TLSSolver::reinitLevelSetField() {
post_pro.exportOnSpace("phi", step, eval_phi, integ_rule, begin, end);
xEvalGradField<xNorm> eval_norm_grad_phi(field);
post_pro.exportOnSpace("norm_grad_phi", step, eval_norm_grad_phi, integ_rule, begin, end);
// post_pro.saveMesh("comp_for_mesh_generation", step, geom.getMesh());
// post_pro.saveField("phi_for_mesh_generation", step, field, geom.getRegion("tls"));
}
void TLSSolver::correctLevelSetField() {
......@@ -295,9 +268,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"));
post_pro.exportOnSpace("transfered_before_phi", Observer::tell<int>("step"), eval_phi, geom.getIntegRuleBasic(0), geom.begin("tls"), geom.end("tls"));
FastMarchingReinit(field,
geom.getRegion("duplicated"),
geom.begin(0, "bnd_fds"), geom.end(0, "bnd_fds"),
geom.begin(0, "int_duplicated"), geom.end(0, "int_duplicated"),
[](const mEntity& v) { return 1.; }, false);
// TransferField(field, field,
// [](xValue<double>* v1, xValue<double>* v2) { v2->setVal(v1->getVal()); },
// geom.begin(0, "duplicated"), geom.end(0, "duplicated"));
deleteField("int_fds");
declareState();
post_pro.exportOnSpace("transfered_phi", Observer::tell<int>("step"), eval_phi, geom.getIntegRuleBasic(0), geom.begin("tls"), geom.end("tls"));
......@@ -325,6 +304,7 @@ void TLSSolver::delocalizeMeanField() {
nonlocality_graph=nullptr;
}
computeModesFastMarching();
geom.deleteBndAndInt();
}
void TLSSolver::computeModesFastMarching() {
......@@ -402,6 +382,15 @@ void TLSSolver::saveStep(int step) {
post_pro.saveField("phi", step, field, geom.getRegion("tls"));
}
void TLSSolver::restoreStep(PreProcessing& pre_pro) {
OldAndCurrent_c::old();
pre_pro.loadField("old_phi", geom.getMesh(), field);
OldAndCurrent_c::current();
pre_pro.loadField("phi", geom.getMesh(), field);
xEvalGradField<xIdentity<xVector> > eval_grad_phi(field);
geom.updateBndAndInt(eval_grad_phi);
}
void TLSSolver::exportStep(int step) {
auto integ_rule=geom.getIntegRulePartition(0);
auto begin=geom.begin("tls");
......@@ -435,8 +424,3 @@ const xEval<double>& TLSSolver::getPhiEval() const {
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"));
// }
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