Commit cc8a097a authored by Benoit LE's avatar Benoit LE

Modifications to take into account templatization of xTensor

parent 9d0cd807
......@@ -125,11 +125,11 @@ void PostProcessing::exportOnTime(const std::string name, const double time, con
sensor_manager.measureAll(eval, time, name);
}
void PostProcessing::exportOnTime(const std::string name, const double time, const xEval<xtensor::xVector>& eval) {
void PostProcessing::exportOnTime(const std::string name, const double time, const xEval<xtensor::xVector<>>& eval) {
sensor_manager.measureAll(eval, time, name);
}
void PostProcessing::exportOnTime(const std::string name, const double time, const xEval<xtensor::xTensor2>& eval) {
void PostProcessing::exportOnTime(const std::string name, const double time, const xEval<xtensor::xTensor2<>>& eval) {
sensor_manager.measureAll(eval, time, name);
}
......
......@@ -77,7 +77,7 @@ public:
}
template<typename ITER>
void exportOnSpace(const std::string export_name, const int step, const xfem::xEval<xtensor::xVector>& eval,
void exportOnSpace(const std::string export_name, const int step, const xfem::xEval<xtensor::xVector<>>& eval,
const xfem::xIntegrationRule& integration_rule, ITER begin, ITER end,
xfem::xEntityToEntity integ2appro = xtool::xIdentity < AOMD::mEntity * >()) {
if(export_manager.toExport(export_name, step, "")) {
......@@ -85,15 +85,15 @@ public:
}
for(int i=0; i<3; ++i) {
if(export_manager.toExport(export_name+"_"+std::to_string(i+1), step, "")) {
xtensor::xExtractCompVector extract(i);
xfem::xEvalUnary<xtensor::xExtractCompVector> evalex(extract, eval);
xtensor::xExtractCompVector<> extract(i);
xfem::xEvalUnary<xtensor::xExtractCompVector<>> evalex(extract, eval);
exportOnSpace(export_name, evalex, integration_rule, begin, end, integ2appro);
}
}
}
template<typename ITER>
void exportOnSpace(const std::string export_name, const int step, const xfem::xEval<xtensor::xTensor2>& eval,
void exportOnSpace(const std::string export_name, const int step, const xfem::xEval<xtensor::xTensor2<>>& eval,
const xfem::xIntegrationRule& integration_rule, ITER begin, ITER end,
xfem::xEntityToEntity integ2appro = xtool::xIdentity < AOMD::mEntity * >()) {
if(export_manager.toExport(export_name, step, "")) {
......@@ -102,18 +102,18 @@ public:
for(int i=0; i<3; ++i) {
for(int j=0; j<3; ++j) {
if(export_manager.toExport(export_name+"_"+std::to_string(i+1)+std::to_string(j+1), step, "")) {
xtensor::xExtractCompTensor extract(i,j);
xfem::xEvalUnary<xtensor::xExtractCompTensor> evalex(extract, eval);
xtensor::xExtractCompTensor<> extract(i,j);
xfem::xEvalUnary<xtensor::xExtractCompTensor<>> evalex(extract, eval);
exportOnSpace(export_name, evalex, integration_rule, begin, end, integ2appro);
}
}
}
if(export_manager.toExport(export_name+"_tr", step, "")) {
xfem::xEvalUnary<xtensor::xTrace> evaltr(eval);
xfem::xEvalUnary<xtensor::xTrace<>> evaltr(eval);
exportOnSpace(export_name, evaltr, integration_rule, begin, end, integ2appro);
}
if(export_manager.toExport(export_name+"_dev", step, "")) {
xfem::xEvalUnary<xtensor::xDeviatoric> evaldev(eval);
xfem::xEvalUnary<xtensor::xDeviatoric<>> evaldev(eval);
exportOnSpace(export_name, evaldev, integration_rule, begin, end, integ2appro);
}
}
......@@ -134,8 +134,8 @@ public:
void exportOnTime(const std::string name, const double time, const double val);
void exportOnTime(const std::string name, const double time, const xfem::xEval<double>& eval);
void exportOnTime(const std::string name, const double time, const xfem::xEval<xtensor::xVector>& eval);
void exportOnTime(const std::string name, const double time, const xfem::xEval<xtensor::xTensor2>& eval);
void exportOnTime(const std::string name, const double time, const xfem::xEval<xtensor::xVector<>>& eval);
void exportOnTime(const std::string name, const double time, const xfem::xEval<xtensor::xTensor2<>>& eval);
void saveMesh(std::string export_name, int step, xfem::xMesh& mesh);
void saveField(std::string export_name, int step, const xfem::xField& field, xfem::xRegion reg);
......
......@@ -93,7 +93,7 @@ void Formulation::transferDispField(double coeff) {
deleteDispField("int_fds");
applyEssentialEnv(coeff);
declareDofs();
xEvalField<xtool::xIdentity<xtensor::xVector> > eval_disp_field(disp_field);
xEvalField<xtool::xIdentity<xtensor::xVector<>> > eval_disp_field(disp_field);
post_pro.exportOnSpace("transfered_disp", Observer::tell<int>("step"), eval_disp_field, geom.getIntegRuleBasic(0), geom.begin(), geom.end());
}
......@@ -177,7 +177,7 @@ void Formulation::addDisp(const xlinalg::xCSRVector& vec) {
void Formulation::writeStrain(std::string name, std::string domain_name) {
OldAndCurrent_c::current();
SetMaterialVariablesVisitor_c<xEvalGradField<xtensor::xSymmetrize>> set_visitor(name, disp_field);
SetMaterialVariablesVisitor_c<xEvalGradField<xtensor::xSymmetrize<>>> set_visitor(name, disp_field);
UpdateMaterialVariablesVisitor_c update_visitor(name);
xeCompositeMaterialVariablesVisitor_c composite_visitor(set_visitor, update_visitor);
VisitMaterialVariablesCommand_c visit_command(composite_visitor, variab_manager);
......@@ -189,20 +189,20 @@ void Formulation::assembleGraphMatrix(LinearSystem& system) {
}
void Formulation::assembleStiffnessMatrix(LinearSystem& system, std::string matrix_name, std::string domain_name) {
NonUniformMaterialSensitivity_c<xtensor::xTensor4> eval_stiffness(matrix_name, variab_manager);
xFormBilinearWithLaw<xGradOperator<xtensor::xSymmetrize>, xEval<xtensor::xTensor4>, xGradOperator<xtensor::xSymmetrize>> form_bilinear(eval_stiffness);
NonUniformMaterialSensitivity_c<xtensor::xTensor4<>> eval_stiffness(matrix_name, variab_manager);
xFormBilinearWithLaw<xGradOperator<xtensor::xSymmetrize<>>, xEval<xtensor::xTensor4<>>, xGradOperator<xtensor::xSymmetrize<>>> form_bilinear(eval_stiffness);
Assemble(form_bilinear, system.getAssembler(), geom.getIntegRuleSmart(), disp_field, disp_field, geom.begin(domain_name), geom.end(domain_name));
}
void Formulation::assembleMassMatrix(LinearSystem& system) {
xUniformMaterialSensitivity<double> eval_density("density");
xFormBilinearWithLaw<xValOperator<xtool::xIdentity<xtensor::xVector>>, xEval<double>, xValOperator<xtool::xIdentity<xtensor::xVector>>> form_bilinear(eval_density);
xFormBilinearWithLaw<xValOperator<xtool::xIdentity<xtensor::xVector<>>>, xEval<double>, xValOperator<xtool::xIdentity<xtensor::xVector<>>>> form_bilinear(eval_density);
Assemble(form_bilinear, system.getAssemblerLumpedEqu(), geom.getIntegRuleSmart(), disp_field, disp_field, geom.begin(), geom.end());
}
void Formulation::assembleInternalForce(LinearSystem& system, std::string vector_name, std::string domain_name) {
NonUniformMaterialSensitivity_c<xtensor::xTensor2> eval_stress(vector_name, variab_manager);
xFormLinearWithLoad<xGradOperator<xtensor::xSymmetrize>, xEval<xtensor::xTensor2> > form_linear(eval_stress);
NonUniformMaterialSensitivity_c<xtensor::xTensor2<>> eval_stress(vector_name, variab_manager);
xFormLinearWithLoad<xGradOperator<xtensor::xSymmetrize<>>, xEval<xtensor::xTensor2<>> > form_linear(eval_stress);
Assemble(form_linear, system.getAssembler(), geom.getIntegRuleSmart(), disp_field, geom.begin(domain_name), geom.end(domain_name));
}
......@@ -215,27 +215,27 @@ void Formulation::assembleNaturalEnv(LinearSystem& system) {
std::string phys=env.Phys;
if(strncmp(phys.c_str(), "TRACTION", 8)==0) {
assert(env.Type==FIX || env.Type==FIX_AND_MEASURE);
xtensor::xVector val;
xtensor::xVector<> val;
if(phys[9]=='X') { val(0)=env.getValue(); }
else if(phys[9]=='Y') { val(1)=env.getValue(); }
else if(phys[9]=='Z') { val(2)=env.getValue(); }
xEvalConstant<xtensor::xVector> flux(val);
xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector> >, xEvalConstant<xtensor::xVector> > lin(flux);
xEvalConstant<xtensor::xVector<>> flux(val);
xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector<>> >, xEvalConstant<xtensor::xVector<>> > lin(flux);
xClassRegion bc(&mesh, env.Entity, env.getDimension());
Assemble(lin, system.getAssembler(), geom.getIntegRuleBasic(bc_integ_order), disp_field, bc.begin(), bc.end(), xUpperAdjacency());
}
else if(phys=="PRESSURE") {
assert(env.Type==FIX || env.Type==FIX_AND_MEASURE);
xEvalNormal eval_normal;
xtool::xScale<xtensor::xVector> scale(env.getValue());
xEvalUnary<xtool::xScale<xtensor::xVector> > eval_pressure(scale, eval_normal);
xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector> >, xEval<xtensor::xVector> > lin(eval_pressure);
xtool::xScale<xtensor::xVector<>> scale(env.getValue());
xEvalUnary<xtool::xScale<xtensor::xVector<>> > eval_pressure(scale, eval_normal);
xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector<>> >, xEval<xtensor::xVector<>> > lin(eval_pressure);
xClassRegion bc(&mesh, env.Entity, env.getDimension());
Assemble(lin, system.getAssembler(), geom.getIntegRuleBasic(bc_integ_order), disp_field, bc.begin(), bc.end(), xUpperAdjacency());
}
else if(strncmp(phys.c_str(), "STRESS", 6)==0) {
assert(env.Type==FIX || env.Type==FIX_AND_MEASURE);
xtensor::xTensor2 val;
xtensor::xTensor2<> val;
if(phys[7]=='X' && phys[8]=='X') { val(0,0)=env.getValue(); }
else if(phys[7]=='X' && phys[8]=='Y') { val(0,1)=env.getValue(); val(1,0)=env.getValue(); }
else if(phys[7]=='X' && phys[8]=='Z') { val(0,2)=env.getValue(); val(2,0)=env.getValue(); }
......@@ -243,9 +243,9 @@ void Formulation::assembleNaturalEnv(LinearSystem& system) {
else if(phys[7]=='Y' && phys[8]=='Z') { val(1,2)=env.getValue(); val(2,1)=env.getValue(); }
else if(phys[7]=='Z' && phys[8]=='Z') { val(2,2)=env.getValue(); }
xEvalNormal eval_normal;
xEvalConstant<xtensor::xTensor2> eval_stress(val);
xEvalBinary<xtool::xMult<xtensor::xTensor2, xtensor::xVector, xtensor::xVector> > eval_traction(eval_stress, eval_normal);
xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector> >, xEval<xtensor::xVector> > lin(eval_traction);
xEvalConstant<xtensor::xTensor2<>> eval_stress(val);
xEvalBinary<xtool::xMult<xtensor::xTensor2<>, xtensor::xVector<>, xtensor::xVector<>> > eval_traction(eval_stress, eval_normal);
xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector<>> >, xEval<xtensor::xVector<>> > lin(eval_traction);
xClassRegion bc(&mesh, env.Entity, env.getDimension());
Assemble(lin, system.getAssembler(), geom.getIntegRuleBasic(bc_integ_order), disp_field, bc.begin(), bc.end(), xUpperAdjacency());
}
......@@ -298,18 +298,18 @@ void Formulation::exportStep(int step) {
post_pro.exportOnSpace(it->first, step, eval, integ_rule, begin, end);
}
for(xTensorsSignature::const_iterator it=tensor_sig.begin_tensor2(); it!=tensor_sig.end_tensor2(); ++it) {
SmoothMaterialVariable_c<xtensor::xTensor2> eval(it->first, variab_manager);
SmoothMaterialVariable_c<xtensor::xTensor2<>> eval(it->first, variab_manager);
post_pro.exportOnSpace(it->first, step, eval, integ_rule, begin, end);
}
for(auto str: {"strain", "stress"}) {
// SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2> eval(str, variab_manager);
// SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2<>> eval(str, variab_manager);
// post_pro.exportOnSpace(str, step, eval, integ_rule, begin, end);
SmoothNonUniformMaterialSensitivity_c<xtensor::xTensor2> eval(str, variab_manager);
SmoothNonUniformMaterialSensitivity_c<xtensor::xTensor2<>> eval(str, variab_manager);
post_pro.exportOnSpace(str, step, eval, geom.getIntegRuleSmart(), begin, end);
}
CustomValueDoubleSelector::which_val=size_disp_value;
xEvalField<xtool::xIdentity<xtensor::xVector> > eval_disp(disp_field);
xEvalField<xtool::xIdentity<xtensor::xVector<>> > eval_disp(disp_field);
post_pro.exportOnSpace("disp", step, eval_disp, integ_rule, begin, end);
post_pro.exportOnSpace("disp_basic", step, eval_disp, geom.getIntegRuleBasic(0), begin, end);
post_pro.exportOnTime("disp", step, eval_disp);
......
......@@ -147,9 +147,9 @@ void DynFormulation::exportStep(int step, double time) {
post_pro.exportOnSpace(str, step, eval, geom.getIntegRuleSmart(), begin, end);
}
for(auto str: {"stress"}) {
// SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2> eval(str, variab_manager);
// SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2<>> eval(str, variab_manager);
// post_pro.exportOnSpace(str, step, eval, integ_rule, begin, end);
SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2> eval(std::string(str)+"_ref", variab_manager);
SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2<>> eval(std::string(str)+"_ref", variab_manager);
post_pro.exportOnSpace(str, step, eval, geom.getIntegRuleSmart(), begin, end);
}
......
......@@ -324,16 +324,16 @@ void QSFormulation::writeCurrentToOld() {
tls_solver.writeCurrentToOldLevelSetField();
}
xtensor::xVector QSFormulation::computeRefForce() {
xtensor::xVector<> QSFormulation::computeRefForce() {
const int nb_dofs = double_manager_meas.size(dofs);
xlinalg::xCSRVector Fint(nb_dofs);
xAssemblerBasic<> assembler_Fint(Fint);
NonUniformMaterialSensitivity_c<xtensor::xTensor2> eval_stress("stress_ref", variab_manager);
xFormLinearWithLoad<xGradOperator<xtensor::xSymmetrize>, xEval<xtensor::xTensor2>> form_linear(eval_stress);
NonUniformMaterialSensitivity_c<xtensor::xTensor2<>> eval_stress("stress_ref", variab_manager);
xFormLinearWithLoad<xGradOperator<xtensor::xSymmetrize<>>, xEval<xtensor::xTensor2<>>> form_linear(eval_stress);
Assemble(form_linear, assembler_Fint, geom.getIntegRuleSmart(), disp_field_meas, geom.begin("bnd_meas"), geom.end("bnd_meas"));
char dir[3]={'X', 'Y', 'Z'};
xtensor::xVector ref_force;
xtensor::xVector<> ref_force;
std::string lis_phys[3]={"DISPLACEMENT_X", "DISPLACEMENT_Y", "DISPLACEMENT_Z"};
xMesh& mesh=geom.getMesh();
......@@ -381,9 +381,9 @@ xtensor::xVector QSFormulation::computeRefForce() {
return ref_force;
}
xtensor::xVector QSFormulation::computeRefDisp() {
xtensor::xVector<> QSFormulation::computeRefDisp() {
xMesh& mesh=geom.getMesh();
xtensor::xVector ref_disp;
xtensor::xVector<> ref_disp;
OldAndCurrent_c::current();
for(xPhysicalEnv::const_iterator it=data.PhysicalEnv->begin(); it!=data.PhysicalEnv->end(); ++it) {
const xEnv& env=*it;
......@@ -397,11 +397,11 @@ xtensor::xVector QSFormulation::computeRefDisp() {
if(phys[9]=='Y') { i=1; }
if(phys[9]=='Z') { i=2; }
xClassRegion bc(&mesh, env.Entity, env.getDimension());
xtensor::xVector vec;
xtensor::xVector<> vec;
vec(i)=1.;
xEvalField<xtool::xIdentity<xtensor::xVector> > eval_disp(disp_field);
xEvalConstant<xtensor::xVector> eval_dir(vec);
xEvalBinary<xtool::xMult<xtensor::xVector, xtensor::xVector, double> > eval(eval_disp, eval_dir);
xEvalField<xtool::xIdentity<xtensor::xVector<>> > eval_disp(disp_field);
xEvalConstant<xtensor::xVector<>> eval_dir(vec);
xEvalBinary<xtool::xMult<xtensor::xVector<>, xtensor::xVector<>, double> > eval(eval_disp, eval_dir);
xIntegrateEvalCommand<xEval<double> > command(eval, ref_disp(i));
ApplyCommandOnIntegrationRule(command, geom.getIntegRuleBasic(), bc.begin(), bc.end());
}
......@@ -433,12 +433,12 @@ bool QSFormulation::exportStep(int step) {
post_pro.exportOnSpace(str, step, eval, geom.getIntegRuleSmart(), begin, end);
}
for(auto str: {"stress_ref"}) {
// SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2> eval(str, variab_manager);
// SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2<>> eval(str, variab_manager);
// post_pro.exportOnSpace(str, step, eval, integ_rule, begin, end);
SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2> eval(str, variab_manager);
SmoothApproNonUniformMaterialSensitivity_c<xtensor::xTensor2<>> eval(str, variab_manager);
post_pro.exportOnSpace(str, step, eval, geom.getIntegRuleSmart(), begin, end);
}
xEvalField<xtool::xIdentity<xtensor::xVector> > eval_disp(disp_field);
xEvalField<xtool::xIdentity<xtensor::xVector<>> > eval_disp(disp_field);
post_pro.exportOnSpace("disp_ref", step, eval_disp, integ_rule, begin, end);
post_pro.exportOnSpace("disp_ref_basic", step, eval_disp, geom.getIntegRuleBasic(0), begin, end);
......@@ -507,7 +507,7 @@ void QSFormulation::saveStep(int step) {
void QSFormulation::restoreStep(PreProcessing& pre_pro) {
pre_pro.loadField("disp", geom.getMesh(), disp_field);
xEvalField<xtool::xIdentity<xtensor::xVector> > eval_disp(disp_field);
xEvalField<xtool::xIdentity<xtensor::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());
......
......@@ -70,8 +70,8 @@ protected:
void clearDispMeasField();
void declareDispMeasField();
xtensor::xVector computeRefForce();
xtensor::xVector computeRefDisp();
xtensor::xVector<> computeRefForce();
xtensor::xVector<> computeRefDisp();
// Returns true if it is the peak step in force/displacement curve.
bool isPeakForce();
......@@ -87,7 +87,7 @@ protected:
xfem::xField disp_field_meas;
double old_load_factor, load_factor, load_factor_variation, max_load_factor, max_force_norm, dissipated_energy;
xtensor::xVector old_ref_force, ref_force, old_ref_disp, ref_disp;
xtensor::xVector<> old_ref_force, ref_force, old_ref_disp, ref_disp;
const double delta_phi_max;
bool is_forced, is_nonlinear, do_anticipated_exit;
};
......
......@@ -414,7 +414,7 @@ void QSFormulRemeshAniso::projectPhiOnMesh(xMesh *mesh_new, std::map<mEntity *,
mEntity *elt = e->get(2,0);
xfem::xGeomElem geo(elt);
Trellis_Util::mPoint CGD = geo.getCDGxyz();
xtensor::xVector toCDG(p,CGD);
xtensor::xVector<> toCDG(p,CGD);
//toCDG.norm();
double epsi = 1.e-3;
......@@ -497,15 +497,15 @@ void QSFormulRemeshAniso::remesh() {
//xMesh* meshForRemeshing=&geom.getMesh();
//Misc evaluators
xEvalField<xtool::xIdentity<xtensor::xVector> > eval_disp(disp_field);
xEvalGradField<xtensor::xSymmetrize> eval_strain(disp_field);
xEvalField<xtool::xIdentity<xtensor::xVector<>> > eval_disp(disp_field);
xEvalGradField<xtensor::xSymmetrize<>> eval_strain(disp_field);
xEvalLevelSet<xtool::xIdentity<double> > eval_lsetAdapt(lsetAdapt);
//Adapters for the evaluators that are defined on geom.getMesh, but exported on meshForRemeshing (backupMesh)
xEvalOnOtherMesh<double> eval_phi_on_backupMesh(tls_solver.getPhiEval(),geom.getMesh() );
xEvalOnOtherMesh<double> eval_damage_on_backupMesh(tls_solver.getDamageEval(),geom.getMesh() );
xEvalOnOtherMesh<xtensor::xVector> eval_disp_on_backupMesh(eval_disp,geom.getMesh() );
xEvalOnOtherMesh<xtensor::xTensor2> eval_strain_on_backupMesh(eval_strain,geom.getMesh() );
xEvalOnOtherMesh<xtensor::xVector<>> eval_disp_on_backupMesh(eval_disp,geom.getMesh() );
xEvalOnOtherMesh<xtensor::xTensor2<>> eval_strain_on_backupMesh(eval_strain,geom.getMesh() );
xEvalOnOtherMesh<double> eval_lsadapt_on_backupMesh(eval_lsetAdapt,geom.getMesh() );
......@@ -519,7 +519,7 @@ void QSFormulRemeshAniso::remesh() {
//postProVTU.exportLevelSet(eval_lsadapt_on_backupMesh,"lsadapt");
postProVTU.exportScalarAtNode(eval_lsadapt_on_backupMesh,"lsadapt");
postProVTU.exportScalarAtNode(eval_damage_on_backupMesh,"Endomagement");
//SmoothNonUniformMaterialSensitivity_c<xtensor::xTensor2> eval("strain", variab_manager);
//SmoothNonUniformMaterialSensitivity_c<xtensor::xTensor2<>> eval("strain", variab_manager);
//TODO
// postProVTU.exportVonMisesAtElem(eval_strain_on_backupMesh,"deformation");
......@@ -612,7 +612,7 @@ void QSFormulRemeshAniso::remesh() {
xReadFromParaview(entree,"deplacement",disp_field,*mesh);
disp_field.getDoubleManager()->PrintForDebug("aft_read_disp.dbg");
xEvalField<xtool::xIdentity<xtensor::xVector> > eval_disp1(disp_field);
xEvalField<xtool::xIdentity<xtensor::xVector<>> > eval_disp1(disp_field);
//xEvalConstant<double> eval_disp_(1.);
post_pro.exportOnSpace("read_disp", Observer::tell<int>("step"), eval_disp1, geom.getIntegRuleSmart(), geom.begin(), geom.end());
......@@ -621,7 +621,7 @@ void QSFormulRemeshAniso::remesh() {
// xexport::xPostProVTU postProVTU2(mesh);
// postProVTU2.openFile("DispRead.vtu");
// xEvalField<xtool::xIdentity<xtensor::xVector> > eval_disp2(disp_field);
// xEvalField<xtool::xIdentity<xtensor::xVector<>> > eval_disp2(disp_field);
// postProVTU2.exportVectorAtNode(eval_disp2,"deplacement");
// postProVTU2.closeFile();
......@@ -688,7 +688,7 @@ void QSFormulRemeshAniso::remesh() {
// tls_solver.updateDomains();
tls_solver.delocalizeMeanField();
xEvalGradField<xtensor::xNorm> eval_norm_grad_phi_after(tls_solver.getField());
xEvalGradField<xtensor::xNorm<>> eval_norm_grad_phi_after(tls_solver.getField());
post_pro.exportOnSpace("norm_grad_phi_after", Observer::tell<int>("step"), eval_norm_grad_phi_after, geom.getIntegRuleBasic(0), geom.begin(), geom.end());
// cout<<"delocalized list size"<<tls_solver.delocalized_list.size()<<endl;
......
......@@ -20,8 +20,8 @@ using namespace xfem;
QSFormulationSuperimposed::QSFormulationSuperimposed(TLSGeom& geom, TLSSolver& tls_solver,
const xData& data, const xParseData& parse_data,
PreProcessing& pre_pro, PostProcessing& post_pro,
const xEval<xtensor::xTensor2>& eval_strain_superimposed,
const xEval<xtensor::xTensor2>& eval_stress_superimposed) :
const xEval<xtensor::xTensor2<>>& eval_strain_superimposed,
const xEval<xtensor::xTensor2<>>& eval_stress_superimposed) :
QSFormulation(geom, tls_solver, data, parse_data, pre_pro, post_pro),
eval_strain_superimposed(eval_strain_superimposed),
eval_stress_superimposed(eval_stress_superimposed)
......@@ -29,9 +29,9 @@ QSFormulationSuperimposed::QSFormulationSuperimposed(TLSGeom& geom, TLSSolver& t
QSFormulationSuperimposed::~QSFormulationSuperimposed() {}
class MyNormal : public xEval<xtensor::xVector> {
class MyNormal : public xEval<xtensor::xVector<>> {
public:
void operator()(const xfem::xGeomElem* appro, const xfem::xGeomElem* integ, xtensor::xVector& vec) const {
void operator()(const xfem::xGeomElem* appro, const xfem::xGeomElem* integ, xtensor::xVector<>& vec) const {
auto pt=integ->getXYZ();
double theta=atan2(pt(1), pt(0));
vec(0)=cos(theta);
......@@ -41,25 +41,25 @@ public:
};
void QSFormulationSuperimposed::assembleNaturalEnv(LinearSystem& system) {
NonUniformMaterialSensitivity_c<xtensor::xTensor4> eval_stiffness("secant_matrix", variab_manager);
xEvalBinary<xtool::xMult<xtensor::xTensor4, xtensor::xTensor2, xtensor::xTensor2>> eval_stress(eval_stiffness, eval_strain_superimposed);
xFormLinearWithLoad<xGradOperator<xtensor::xSymmetrize>, xEval<xtensor::xTensor2> > form_linear(eval_stress);
NonUniformMaterialSensitivity_c<xtensor::xTensor4<>> eval_stiffness("secant_matrix", variab_manager);
xEvalBinary<xtool::xMult<xtensor::xTensor4<>, xtensor::xTensor2<>, xtensor::xTensor2<>>> eval_stress(eval_stiffness, eval_strain_superimposed);
xFormLinearWithLoad<xGradOperator<xtensor::xSymmetrize<>>, xEval<xtensor::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());
system.setAssemblerCoeff(1.);
MyNormal eval_normal;
xEvalBinary<xtool::xMult<xtensor::xTensor2, xtensor::xVector, xtensor::xVector> > eval_traction(eval_stress_superimposed, eval_normal);
xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector> >, xEval<xtensor::xVector> > lin(eval_traction);
xEvalBinary<xtool::xMult<xtensor::xTensor2<>, xtensor::xVector<>, xtensor::xVector<>> > eval_traction(eval_stress_superimposed, eval_normal);
xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector<>> >, xEval<xtensor::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());
}
void QSFormulationSuperimposed::writeStrain(std::string name, std::string domain_name) {
OldAndCurrent_c::current();
xEvalGradField<xtensor::xSymmetrize> eval_strain_tilde(disp_field);
xEvalBinary<std::plus<xtensor::xTensor2> > eval_strain(eval_strain_tilde, eval_strain_superimposed);
SetMaterialVariablesVisitor_c<xEvalBinary<std::plus<xtensor::xTensor2>>> set_visitor(name, eval_strain);
xEvalGradField<xtensor::xSymmetrize<>> eval_strain_tilde(disp_field);
xEvalBinary<std::plus<xtensor::xTensor2<>> > eval_strain(eval_strain_tilde, eval_strain_superimposed);
SetMaterialVariablesVisitor_c<xEvalBinary<std::plus<xtensor::xTensor2<>>>> set_visitor(name, eval_strain);
UpdateMaterialVariablesVisitor_c update_visitor(name);
xeCompositeMaterialVariablesVisitor_c composite_visitor(set_visitor, update_visitor);
VisitMaterialVariablesCommand_c visit_command(composite_visitor, variab_manager);
......
......@@ -19,7 +19,7 @@ public:
QSFormulationSuperimposed(TLSGeom&, TLSSolver&,
const xfem::xData&, const xParseData&,
PreProcessing&, PostProcessing&,
const xfem::xEval<xtensor::xTensor2>&, const xfem::xEval<xtensor::xTensor2>&);
const xfem::xEval<xtensor::xTensor2<>>&, const xfem::xEval<xtensor::xTensor2<>>&);
virtual ~QSFormulationSuperimposed();
// Assembles external force vector by reading main.dat file and
......@@ -29,8 +29,8 @@ public:
// Writes strain given by std::string using disp field symmetric gradient.
virtual void writeStrain(std::string, std::string="all");
private:
const xfem::xEval<xtensor::xTensor2>& eval_strain_superimposed;
const xfem::xEval<xtensor::xTensor2>& eval_stress_superimposed;
const xfem::xEval<xtensor::xTensor2<>>& eval_strain_superimposed;
const xfem::xEval<xtensor::xTensor2<>>& eval_stress_superimposed;
};
#endif
......@@ -100,7 +100,7 @@ void ReadData(std::string main_file_name, xfem::xData& data, xParseData& parse_d
parse_data.registerListString("export_manager", default_export_list);
default_export_list.clear();
parse_data.registerListString("export_sensors_label", default_export_list);
std::list<xtensor::xVector> default_export_point;
std::list<xtensor::xVector<>> default_export_point;
parse_data.registerListVector("export_sensors_point", default_export_point);
data.ReadInfo(main_file_name.c_str());
......
......@@ -38,10 +38,10 @@ void ElastoDam::checkProperties() {
}
}
void ElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTensor2& sensitivity) {
void ElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTensor2<>& sensitivity) {
if (phys_token == "stress_ref") {
const double& d=curr->scalar("d");
const xtensor::xTensor2& eps_ref=curr->tensor2("eps_ref");
const xtensor::xTensor2<>& eps_ref=curr->tensor2("eps_ref");
double tr_ref=eps_ref.trace();
sensitivity=eps_ref*(2.*mu*(1.-d))+identity*(tr_ref*lambda*(1.-d));
}
......@@ -64,7 +64,7 @@ void ElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTensor2&
void ElastoDam::sensitivityTo(const std::string& phys_token, double& sensitivity) {
if (phys_token == "a") {
const xtensor::xTensor2& eps_ref=curr->tensor2("eps_ref");
const xtensor::xTensor2<>& eps_ref=curr->tensor2("eps_ref");
double tr_ref=eps_ref.trace();
sensitivity=mu*eps_ref.contract(eps_ref)+0.5*lambda*tr_ref*tr_ref;
}
......@@ -78,8 +78,8 @@ void ElastoDam::sensitivityTo(const std::string& phys_token, double& sensitivity
sensitivity=Y_c*hardening_function->getGrad(curr->scalar("d"));
}
else if (phys_token == "beta") {
const xtensor::xTensor2& eps_ref=curr->tensor2("eps_ref");
xtensor::xTensor2 tmp=eps_ref*(2.*mu)+identity*(eps_ref.trace()*lambda);
const xtensor::xTensor2<>& eps_ref=curr->tensor2("eps_ref");
xtensor::xTensor2<> tmp=eps_ref*(2.*mu)+identity*(eps_ref.trace()*lambda);
double load_factor;
Observer::tell("load_factor", load_factor);
sensitivity=tmp.contract(eps_ref)*load_factor;
......@@ -117,7 +117,7 @@ void ElastoDam::sensitivityTo(const std::string& phys_token, double& sensitivity
}
}
void ElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTensor4& sensitivity) {
void ElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTensor4<>& sensitivity) {
if(phys_token == "elastic_matrix") {
sensitivity = elastic_matrix;
}
......@@ -139,7 +139,7 @@ void ElastoDam::computeCurrentState() {
void ElastoDam::computeCurrentState(std::string name) {
if(name == "eps_ref") {
if(plane_state==PLANE_STRESS) {
xtensor::xTensor2& eps_ref = curr->tensor2(name);
xtensor::xTensor2<>& eps_ref = curr->tensor2(name);
eps_ref(2,2) = (eps_ref(0,0)+eps_ref(1,1))*(nu/(nu-1.));
}
}
......@@ -194,17 +194,17 @@ void AsymElastoDam::checkProperties() {
}
}
void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTensor2& sensitivity) {
void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTensor2<>& sensitivity) {
if (phys_token == "stress_ref") {
const double& d=curr->scalar("d");
const xtensor::xTensor2& eps_ref=curr->tensor2("eps_ref");
const xtensor::xTensor2<>& eps_ref=curr->tensor2("eps_ref");
const double tr_ref=eps_ref.trace();
if(d==0.) {
sensitivity = eps_ref*(2.*mu)+identity*(tr_ref*lambda);
return;
}
const xtensor::xTensor2& eps_ref_pos=curr->tensor2("eps_ref_pos");
const xtensor::xTensor2& eps_ref_neg=curr->tensor2("eps_ref_neg");
const xtensor::xTensor2<>& eps_ref_pos=curr->tensor2("eps_ref_pos");
const xtensor::xTensor2<>& eps_ref_neg=curr->tensor2("eps_ref_neg");
const double tr_ref_pos=std::max(0., tr_ref);
const double tr_ref_neg=std::min(0., tr_ref);
sensitivity = eps_ref_pos*(2.*mu*(1.-d)) + identity*(tr_ref_pos*lambda*(1.-d)) + eps_ref_neg*(2.*mu*(1.-h*d)) + identity*(tr_ref_neg*lambda*(1.-h*d));
......@@ -228,8 +228,8 @@ void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTens
void AsymElastoDam::sensitivityTo(const std::string& phys_token, double& sensitivity) {
if (phys_token == "a") {
const xtensor::xTensor2& eps_ref_pos=curr->tensor2("eps_ref_pos");
const xtensor::xTensor2& eps_ref_neg=curr->tensor2("eps_ref_neg");
const xtensor::xTensor2<>& eps_ref_pos=curr->tensor2("eps_ref_pos");
const xtensor::xTensor2<>& eps_ref_neg=curr->tensor2("eps_ref_neg");
const double tr_ref=curr->tensor2("eps_ref").trace();
const double tr_ref_pos=std::max(0., tr_ref);
const double tr_ref_neg=std::min(0., tr_ref);
......@@ -245,13 +245,13 @@ void AsymElastoDam::sensitivityTo(const std::string& phys_token, double& sensiti
sensitivity=Y_c*hardening_function->getGrad(curr->scalar("d"));
}
else if (phys_token == "beta") {
const xtensor::xTensor2& eps_ref=curr->tensor2("eps_ref");
const xtensor::xTensor2& eps_ref_pos=curr->tensor2("eps_ref_pos");
const xtensor::xTensor2& eps_ref_neg=curr->tensor2("eps_ref_neg");
const xtensor::xTensor2<>& eps_ref=curr->tensor2("eps_ref");
const xtensor::xTensor2<>& eps_ref_pos=curr->tensor2("eps_ref_pos");
const xtensor::xTensor2<>& eps_ref_neg=curr->tensor2("eps_ref_neg");
const double tr_ref=eps_ref.trace();
const double tr_ref_pos=std::max(0., tr_ref);
const double tr_ref_neg=std::min(0., tr_ref);
xtensor::xTensor2 tmp=(eps_ref_pos+eps_ref_neg*h)*(2.*mu)+identity*((tr_ref_pos+tr_ref_neg*h)*lambda);
xtensor::xTensor2<> tmp=(eps_ref_pos+eps_ref_neg*h)*(2.*mu)+identity*((tr_ref_pos+tr_ref_neg*h)*lambda);
double load_factor;
Observer::tell("load_factor", load_factor);
sensitivity=tmp.contract(eps_ref)*load_factor;
......@@ -286,7 +286,7 @@ void AsymElastoDam::sensitivityTo(const std::string& phys_token, double& sensiti
}
}
void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTensor4& sensitivity) {
void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTensor4<>& sensitivity) {
if(phys_token == "elastic_matrix") {
sensitivity=elastic_matrix;
}
......@@ -296,8 +296,8 @@ void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTenso
sensitivity=elastic_matrix;
return;
}
const xtensor::xTensor4& proj_pos_old = curr->tensor4("eps_ref_proj_pos");
const xtensor::xTensor4& delta_proj_pos_old = curr->tensor4("eps_ref_proj_pos_delta");
const xtensor::xTensor4<>& proj_pos_old = curr->tensor4("eps_ref_proj_pos");
const xtensor::xTensor4<>& delta_proj_pos_old = curr->tensor4("eps_ref_proj_pos_delta");
const double& tr_ref=curr->tensor2("eps_ref").trace();
double old_load_factor, load_factor, load_factor_variation;
......@@ -308,7 +308,7 @@ void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTenso
// double ratio=0.;
double ratio=std::fabs(load_factor_variation)/std::fabs(load_factor-old_load_factor)*(1.-curr->scalar("d"));
xtensor::xTensor4 proj_pos_pred=proj_pos_old+delta_proj_pos_old*ratio;
xtensor::xTensor4<> proj_pos_pred=proj_pos_old+delta_proj_pos_old*ratio;
sensitivity=proj_pos_pred*(2.*mu*(1.-d))+(id3-proj_pos_pred)*(2.*mu*(1.-h*d));
if(tr_ref>=0.) {
......@@ -321,7 +321,7 @@ void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTenso
else if(phys_token == "tangent_matrix") {
if(plane_state==PLANE_STRAIN || plane_state==PLANE_STRESS) {
const double& d=curr->scalar("d");
const xtensor::xTensor2& eps_ref=curr->tensor2("eps_ref");
const xtensor::xTensor2<>& eps_ref=curr->tensor2("eps_ref");
// (x,y) plane eigen values computation
double tr=eps_ref(0,0)+eps_ref(1,1);
......@@ -370,14 +370,14 @@ void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTenso
// to understand comments, see TLS first article.
// first term : dedeps^T * secderiv_pot * dedeps
xtensor::xTensor2 potentiel_sderivative;
xtensor::xTensor2<> potentiel_sderivative;
potentiel_sderivative(0,0) = 2.*mu*c1+lambda_c_tilde;
potentiel_sderivative(0,1) = lambda_c_tilde;
potentiel_sderivative(1,1) = 2.*mu*c2+lambda_c_tilde;
potentiel_sderivative(1,0) = potentiel_sderivative(0,1);
// derivative eigenvalues w.r.t. (eps_11, eps_22, 2eps_12)
xtensor::xTensor2 dedeps;
xtensor::xTensor2<> dedeps;
dedeps(0,0) = 1+coss;
dedeps(0,1) = 1-coss;
dedeps(0,2) = sinn;
......@@ -386,14 +386,14 @@ void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTenso
dedeps(1,2) = -sinn;
dedeps*=0.5;
xtensor::xTensor2 dedeps_tr(dedeps);
xtensor::xTensor2<> dedeps_tr(dedeps);
dedeps_tr.transpose();
xtensor::xTensor2 tgt;
xtensor::xTensor2<> tgt;
tgt = dedeps_tr*potentiel_sderivative*dedeps;
// second and third terms
xtensor::xTensor2 M;
xtensor::xTensor2<> M;
double cos2 = coss*coss;
double sin2 = sinn*sinn;
double sc = sinn*coss;
......@@ -413,7 +413,7 @@ void AsymElastoDam::sensitivityTo(const std::string& phys_token, xtensor::xTenso
else
tgt += M*(c1+c2);
// write in xtensor::xTensor4
// write in xtensor::xTensor4<>
for(int i=0; i<3; ++i) {
for(int j=0; j<3; ++j) {
for(int k=0; k<3; ++k) {
......@@ -460,8 +460,8 @@ void AsymElastoDam::computeCurrentState(std::string name) {
if(name == "eps_ref") {
std::string mat_class=properties.astring("MATERIAL_CLASS");
if(mat_class=="elastic_damage_nonlinear") {
xtensor::xTensor2& eps_ref = curr->tensor2("eps_ref");
xtensor::xVector v[3];
xtensor::xTensor2<>& eps_ref = curr->tensor2("eps_ref");
xtensor::xVector<> v[3];
double e[3];
if(plane_state==PLANE_STRAIN) {
eps_ref.eigen2d(v, e);
......@@ -507,30 +507,30 @@ void AsymElastoDam::computeCurrentState(std::string name) {
else {
eps_ref.eigen(v, e);
}
xtensor::xTensor2& eps_ref_pos = curr->tensor2("eps_ref_pos");
xtensor::xTensor2& eps_ref_neg = curr->tensor2("eps_ref_neg");
xtensor::xTensor2<>& eps_ref_pos = curr->tensor2("eps_ref_pos");
xtensor::xTensor2<>& eps_ref_neg = curr->tensor2("eps_ref_neg");
eps_ref_pos=std::max(0.,e[0])*tensor_product(v[0],v[0])+
std::max(0.,e[1])*tensor_product(v[1],v[1])+
std::max(0.,e[2])*tensor_product(v[2],v[2]);
eps_ref_neg=eps_ref-eps_ref_pos;
}
else if(mat_class=="elastic_damage") {
const xtensor::xTensor2& eps_ref = curr->tensor2("eps_ref");
xtensor::xVector v[3];
const xtensor::xTensor2<>& eps_ref = curr->tensor2("eps_ref");
xtensor::xVector<> v[3];
double e[3];
if(plane_state==PLANE_STRAIN) {
eps_ref.eigen2d(v, e);
e[2]=0.;
}
xtensor::xTensor4& proj_pos = curr->tensor4("eps_ref_proj_pos");
xtensor::xTensor2& eps_ref_pos = curr->tensor2("eps_ref_pos");
xtensor::xTensor2& eps_ref_neg = curr->tensor2("eps_ref_neg");
xtensor::xTensor4<>& proj_pos = curr->tensor4("eps_ref_proj_pos");
xtensor::xTensor2<>& eps_ref_pos = curr->tensor2("eps_ref_pos");