Commit 89ca0c12 authored by Alexis SALZMAN's avatar Alexis SALZMAN

Mainly adapting to new FM interface to treat GAMAIN and avoid prb

dam_pull_out
============
Adapting to para version

FastMarchingInterface,TLSGeom,TLSSolver
=======================================
Modify groups (sm, sub mesh) to enter FM with known fixed and volatile
values. This is following new FM interface. It correct bug with previous
interface where progression of the front was stopped by re-initialisation
process (all new variation of phi removed by fast marching).
Now buildBndIntNonlocal (TLSGEOM) is creating 2 set (bnd_part_nonlocal
and bnd_nonlocal) in replacement of the original bnd_nonlocal. Those
sets split boundary of non local zone in :
* boundary with local mesh element (bnd_nonlocal)
* boundary with part boundary (bnd_part_nonlocal)
Has before it is do_bnd_in that control if we create bnd_part_nonlocal or
not by calling GetBoundary or GetInteriorBoundary.
In GetBoundary a dirty hack is used to identify part boundary and remove
Gama_c.
When do_bnd_in is 0 bnd_part_nonlocal exist but is empty.

For FM, bnd_part_nonlocal is used has known volatile value. Says those
part boundary value may or may not be re-initialized by FM. In the
latter case those values are used for velocity computation as new extra
modes with their values unmodified. This imply that on part
boundary the level set live its own life when phi gradient enter the
part.

Unfortunately as can be seen, for now, in new dam_gamain test, setting
bnd_part_nonlocal as volatile values may stop progression of the front
just like before when all modified values are reset by FM.

N. Chevaugeon suggested to create only volatile values for non modified
values (values not updated during current time step). This look like a
promising solution because removing these values from bnd_part_nonlocal
and put them in bnd_nonlocal force them to be still present in velocity
computation. It ensure that variation of phi is not destroyed. But it
does not insure anymore that phi gradient norm is 1 around those values.
This somehow must be the case when gradient direction start to be collinear
to part boundary. In that case epsilon_ratio was supposed to do the job
but for rather blunt variation of phi it is not clear if this will fits
ours needs. So filtering like that bnd_part_nonlocal seems interesting.
TODO

About epsilon_ratio it have been forced to zero because really strange
result where obtain with it (non same result in case where
bnd_part_nonlocal is empty, non deterministic computation !???).
To be analyzed.

dam_gamain
==========
This test case is to check if and how gamin feature works.
It is a simple plate with a shear loading that force creation of crack
with an angle of -45 degres (x axes). A default( removing of 2 elements)
is placed in the diagonal of the plate a little closer to left upper
corner. Mesh is regular and chosen to intentionally offers a line of
nodes from default to corner. lc is chosen rather big so that number of
element per lc is important.

Without gamin non-local zone progress toward corner. When reaching
part boundary modification of phi on it are clearly done by FM re-init only.
It arrive to a point where bluntly in on step all corner is extra damaged (phi>lc).
With gamain when reaching part boundary modification comes clearly from
phi update in far more slower way. Crack continue on the diagonal.
Unfortunately process is stopped by the case of delta-phi removed by FM.
(see above for perspectives)

In both case it is not clear how boundary condition have to be handle in
omega+ zone. Without treatment like here it impose force on almost
damaged material which induce maybe other estimate strain. To be check.
parent 10c05fbd
......@@ -46,15 +46,17 @@ public:
private :
};
template<typename ITERTRIALKNOWN,
typename ITERTOFIND,
typename ITERKNOWN>
template< typename ITERKNOWN
,typename ITERKNOWNV
,typename ITERTOFIND
,typename ITERKNOWNT>
void FastMarchingReinit(std::function<double (AOMD::mVertex*)> get_val,
std::function<void (AOMD::mVertex*, double)> set_val,
const xfem::xRegion& region,
ITERTRIALKNOWN trial_known_begin, ITERTRIALKNOWN trial_known_end,
ITERKNOWN known_begin, ITERKNOWN known_end,
ITERKNOWNV known_volatil_begin, ITERKNOWNV known_volatil_end,
ITERTOFIND tofind_begin, ITERTOFIND tofind_end,
ITERKNOWN known_it,
ITERKNOWNT known_it,
double epsilon_ratio,
std::function<double (const AOMD::mVertex&)> f_func=[](const AOMD::mVertex&){return 1.;}) {
#ifndef HAVE_FASTMARCHING
......@@ -64,17 +66,23 @@ void FastMarchingReinit(std::function<double (AOMD::mVertex*)> get_val,
meshinterfacexRegion mi(region);
entitystorage<meshinterfacexRegion, AOMD::mVertex, double> ls(mi);
entitystorage<meshinterfacexRegion, AOMD::mVertex, vector3d<double> > gls(mi);
for(auto it=trial_known_begin; it!=trial_known_end; ++it) {
for(auto it=known_begin; it!=known_end; ++it) {
auto v=static_cast<AOMD::mVertex*>(*it);
ls.set(*v, get_val(v));
}
for(auto it=known_volatil_begin; it!=known_volatil_end; ++it) {
auto v=static_cast<AOMD::mVertex*>(*it);
ls.set(*v, get_val(v));
}
for(auto it=tofind_begin; it!=tofind_end; ++it) {
ls.set(*static_cast<AOMD::mVertex*>(*it), std::numeric_limits<double>::max());
}
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, epsilon_ratio, gls);
entitytovertexiteratorconvertor<ITERKNOWN> trial_known_conv_begin(known_begin);
entitytovertexiteratorconvertor<ITERKNOWN> trial_known_conv_end(known_end);
entitytovertexiteratorconvertor<ITERKNOWNV> trial_known_volatil_conv_begin(known_volatil_begin);
entitytovertexiteratorconvertor<ITERKNOWNV> trial_known_volatil_conv_end(known_volatil_end);
entitytovertexiteratorconvertor<ITERKNOWNT> known_conv_it(known_it);
fmeik(mi, ls, trial_known_conv_begin, trial_known_conv_end,trial_known_volatil_conv_begin, trial_known_volatil_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;
......@@ -115,12 +123,15 @@ void FastMarchingModeExtension(std::function<double (AOMD::mVertex*)> get_val,
for(auto it=tofind_begin; it!=tofind_end; ++it) {
ls.set(*static_cast<AOMD::mVertex*>(*it), std::numeric_limits<double>::max());
}
const std::set<const AOMD::mVertex*> empty;
entitytovertexiteratorconvertor<ITERKNOWN> known_conv_begin(known_begin);
entitytovertexiteratorconvertor<ITERKNOWN> known_conv_end(known_end);
entitytovertexiteratorconvertor< std::set<const AOMD::mVertex*>::const_iterator> known_volatil_conv_begin(empty.begin());
entitytovertexiteratorconvertor< std::set<const AOMD::mVertex*>::const_iterator> known_volatil_conv_end(empty.end());
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, epsilon_ratio, gls, vn);
fmeik(mi, ls, known_conv_begin, known_conv_end,known_volatil_conv_begin, known_volatil_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);
......
......@@ -59,13 +59,27 @@ void DelAdjacencies(mEntity* e) {
}
}
}
// coment to be changed
// Boundary (sm_bnd) of a part (sm_dom), subset of the closure of the open domain (interface with complement and domain boundary).
void GetBoundary(const xSubMesh& sm_dom, xSubMesh& sm_bnd) {
void GetBoundary(const xSubMesh& sm_dom, xSubMesh& sm_bnd, xSubMesh& sm_bnd_part) {
const int dim=sm_dom.dim();
xAcceptInSubMesh is_in_dom(sm_dom);
for(auto it=sm_dom.begin(dim-1); it!=sm_dom.end(dim-1); ++it) {
auto e=*it;
if(e->size(dim)==1 || !is_in_dom(e->get(dim,0)) || !is_in_dom(e->get(dim,1))) {
if(e->size(dim)==1 )
{
pGEntity g = e->getClassification();
// durty hack ! Don't know where this append but classification for newlly created cut element seem
// to use a classification from face on edges. You can see that with this output :
// cout<<"grrrrrr "<<GEN_tag(g)<< " "<<GEN_type(g)<<" "<<e->getLevel()<<endl;
// We used it to select only true boundary
if (GEN_type(g)==1)
{
sm_bnd_part.add(e);
AddNeighbors(e, sm_bnd_part, 0);
}
}
else if(!is_in_dom(e->get(dim,0)) || !is_in_dom(e->get(dim,1))) {
sm_bnd.add(e);
AddNeighbors(e, sm_bnd, 0);
}
......@@ -93,6 +107,16 @@ void GetComplement(const xSubMesh& sm_dom, const xSubMesh& sm_part, xSubMesh& sm
}
}
}
void GetComplement(const xSubMesh& sm_dom, const xSubMesh& sm_part1, const xSubMesh& sm_part2, xSubMesh& sm_compl, int what) {
xAcceptInSubMesh is_in_part1(sm_part1);
xAcceptInSubMesh is_in_part2(sm_part2);
for(auto it=sm_dom.begin(what); it!=sm_dom.end(what); ++it) {
auto e=*it;
if(!is_in_part1(e) && !is_in_part2(e)) {
sm_compl.add(e);
}
}
}
class DelUpdateAdj {
public:
DelUpdateAdj(xMesh& mesh) :
......@@ -580,19 +604,21 @@ void TLSGeom::buildBndIntNonlocal() {
auto& mesh=getMesh();
auto& sm_nonlocal=mesh.getSubMesh("nonlocal");
auto& sm_bnd_nonlocal=mesh.createSubMesh("bnd_nonlocal");
auto& sm_bnd_part_nonlocal=mesh.createSubMesh("bnd_part_nonlocal");
auto& sm_int_nonlocal=mesh.createSubMesh("int_nonlocal");
if(do_bnd_in) {
GetBoundary(sm_nonlocal, sm_bnd_nonlocal);
GetBoundary(sm_nonlocal, sm_bnd_nonlocal,sm_bnd_part_nonlocal);
}
else {
GetInteriorBoundary(sm_nonlocal, sm_bnd_nonlocal);
}
GetComplement(sm_nonlocal, sm_bnd_nonlocal, sm_int_nonlocal, 0);
GetComplement(sm_nonlocal, sm_bnd_nonlocal, sm_bnd_part_nonlocal, sm_int_nonlocal, 0);
xEvalConstant<double> eval_one(1.);
xIntegrationRuleBasic integ_rule(0);
auto step=Observer::tell<int>("step");
post_pro.exportOnSpace("sm_bnd_nonlocal", step, eval_one, integ_rule, sm_bnd_nonlocal.begin(), sm_bnd_nonlocal.end());
post_pro.exportOnSpace("sm_bnd_part_nonlocal", step, eval_one, integ_rule, sm_bnd_part_nonlocal.begin(), sm_bnd_part_nonlocal.end());
post_pro.exportOnSpace("sm_int_nonlocal", step, eval_one, integ_rule, sm_int_nonlocal.begin(), sm_int_nonlocal.end());
}
......@@ -773,6 +799,7 @@ void TLSGeom::buildLinearAndNonLinear(const xEval<double>& eval_phi) {
void TLSGeom::deleteBndIntNonlocal() {
auto& mesh=getMesh();
mesh.deleteSubMesh("bnd_nonlocal");
mesh.deleteSubMesh("bnd_part_nonlocal");
mesh.deleteSubMesh("int_nonlocal");
}
......
......@@ -59,7 +59,9 @@ TLSSolver::TLSSolver(TLSGeom& geom, const xParseData& parse_data,
xAcceptAll(), xAcceptAll()),
do_crack_cut(static_cast<bool>(parse_data.getInt("do_crack_cut"))),
is_newly_delocalized(false),
epsilon_ratio(static_cast<bool>(parse_data.getInt("do_bnd_in"))?1.e-6:0.) {
//epsilon_ratio(static_cast<bool>(parse_data.getInt("do_bnd_in"))?1.e-6:0.) {
//epsilon_ratio(1.e-6) {
epsilon_ratio(0.) {
space_factory->setSpaceProductOrder(1);
field.insert(space_factory->getSpace("fem"));
......@@ -168,7 +170,8 @@ void TLSSolver::swapOldAndCurrentLevelSetField() {
}
void TLSSolver::updateLevelSetField() {
// exportMeanField("delta_phi", Observer::tell<int>("step")); // time consuming
post_pro.exportOnSpace("phi_before_update", Observer::tell<int>("step"), eval_phi, geom.getIntegRuleBasic(0), geom.begin("tls"), geom.end("tls"));
exportMeanField("delta_phi", Observer::tell<int>("step")); // time consuming
// PARA : nonlocality_matrix will change this call in future
// PARA : see commit tagged WITHMODIF
ConvertMeanVectorModalToNodal(nonlocality_matrix, getMeanVector("delta_phi"), buf_nodal_vec);
......@@ -184,6 +187,7 @@ void TLSSolver::reinitLevelSetField() {
geom.buildBndIntNonlocal();
auto nonlocal=geom.getRegion("nonlocal");
auto bnd_nonlocal=geom.getRegion("bnd_nonlocal");
auto bnd_part_nonlocal=geom.getRegion("bnd_part_nonlocal");
auto int_nonlocal=geom.getRegion("int_nonlocal");
std::vector<mEntity*> bnd_in_nonlocal_vec;
bnd_in_nonlocal_vec.reserve(bnd_nonlocal.size(0));
......@@ -191,6 +195,7 @@ void TLSSolver::reinitLevelSetField() {
[this](mVertex* v, double val){ this->field.setVal(v, val); },
nonlocal,
bnd_nonlocal.begin(0), bnd_nonlocal.end(0),
bnd_part_nonlocal.begin(0), bnd_part_nonlocal.end(0),
int_nonlocal.begin(0), int_nonlocal.end(0),
std::back_inserter(bnd_in_nonlocal_vec), epsilon_ratio);
correctLevelSetField();
......
cmake_minimum_required(VERSION 2.6)
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../../../Util/cmakeUtil/)
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../../../Xfiles/Util/cmakeUtil/)
include(common_functions)
include(devel_common_functions)
......@@ -10,50 +10,55 @@ enable_testing()
# Please, sort test cases so that it is easier to maintain
if(TLSDUCTILE_TEST_DAM)
file(GLOB LIST
list(APPEND LISTTEST
# dam_3holes
# dam_3pt_bending
# dam_crack
# dam_hole
# dam_homogeneous
dam_notch
dam_gamain
# dam_pressured_cylinder
# dam_pull_out
dam_pull_out
)
create_tests_from_list(${LIST})
endif()
if(TLSDUCTILE_TEST_DYN)
file(GLOB LIST
list(APPEND LISTTEST
# dyn_branching
# dyn_edge_on_impact
dyn_kalthoff
# dyn_pressured_cylinder
)
create_tests_from_list(${LIST})
endif()
if(TLSDUCTILE_TEST_PLAST)
file(GLOB LIST
list(APPEND LISTTEST
# plast_flounder_eyes
# plast_full_beam
# plast_quarter_beam
# plast_two_holes
)
create_tests_from_list(${LIST})
endif()
if(TLSDUCTILE_TEST_ANIMESH_DAM)
file(GLOB LIST
list(APPEND LISTTEST
# animesh_dam_notch
)
create_tests_from_list(${LIST})
endif()
if(TLSDUCTILE_TEST_INFINITE_DAM)
file(GLOB LIST
list(APPEND LISTTEST
# infinite_dam_homogeneous
# infinite_dam_notch
)
create_tests_from_list(${LIST})
endif()
set(LIST "")
foreach (v ${LISTTEST} )
list(APPEND LIST
${CMAKE_CURRENT_SOURCE_DIR}/${v}
)
endforeach(v)
create_tests_from_list(${LIST})
enable_testing()
add_test(dam_gamain ${CMAKE_CURRENT_BINARY_DIR}/dam_gamain ${CMAKE_CURRENT_BINARY_DIR}/data/main.dat)
add_test(ndiff_dam_gamain dam_gamain ${DEVROOT}/Util/cmakeUtil/test_ndiff.sh )
set_tests_properties(ndiff_dam_gamain PROPERTIES DEPENDS dam_gamain )
integ_order = 3
disp_space_type = PolyLagrange
disp_space_dim = V2Dxy
disp_space_order = 1
disp_bc_integ_order = 1
mean_field_space_type = PolyLagrange
damage_shape = Poly2Revert
lc = 0.6
delta_phi_max_ratio = 0.99
do_forced_delocalization = 0
do_crack_cut = 1
do_bnd_in = 1
fit_to_vertices_ratio = -1.e-2
nb_step_max = 999
export_manager = {
disp 10
phi 1
phi_after_update 1
phi_before_update 1
sm_local 1
sm_bnd_nonlocal 1
sm_bnd_part_nonlocal 1
sm_int_nonlocal 1
sm_bnd_in_nonlocal 1
sm_int_in_nonlocal 1
sm_nonlocal 1
damage 10 }
export_sensors_label = { load_factor 1 force_y 1 dissipated_energy 1 }
export_sensors_point = { }
MESH_FILE_TYPE = msh
PROCEDURE_PARAM_FILE = data/info.dat
MESH_FILE = data/mesh.msh
ZONE 26 = { MAT_CLASS = ElastoDam MAT_PARAM = data/mate.mat }
#teta = 45
# Sig*(-sin(teta)*{ cos(teta), sin(teta) }
BC_LINE 11 ={ TRACTION_X FIX = -0.5e6
TRACTION_Y FIX = -0.5e6 }
# Sig*(+cos(teta)*{ cos(teta), sin(teta) }
BC_LINE 12 ={ TRACTION_X FIX = 0.5e6
TRACTION_Y FIX = 0.5e6 }
# Sig*(+sin(teta)*{ cos(teta), sin(teta) }
BC_LINE 13 ={ TRACTION_X FIX = 0.5e6
TRACTION_Y FIX = 0.5e6 }
# Sig*(-cos(teta)*{ cos(teta), sin(teta) }
BC_LINE 14 ={ TRACTION_X FIX = -0.5e6
TRACTION_Y FIX = -0.5e6 }
BC_POINT 8 ={ DISPLACEMENT_Y FIX = 0.0}
BC_POINT 7 ={ DISPLACEMENT_X FIX = 0.0
DISPLACEMENT_Y FIX = 0.0 }
MATERIAL_CLASS = elastic_damage
NAME = truc
PLANE_STATE = plane_strain
YOUNG_MODULUS = 210.e9
POISSON_RATIO = 0.3
Y_CRIT = 1.32e5
HARDENING = exponential
HARDENING_COEFF = 4.
This diff is collapsed.
/*
This source code is subject to non-permissive licence,
see the TLSDuctile/LICENSE file for conditions.
*/
#include "TLSDuctile.h"
int main(int argc, char* argv[]) {
eXlibris_tools::xMPIEnv::init(argc,argv);
RegisterMaterials();
std::string data_filename, archive_filename;
ParseArgs(argc, argv, data_filename, archive_filename);
xfem::xData data;
xParseData parse_data;
ReadData(data_filename, data, parse_data);
xfem::xExportGmshAscii pexport;
PreProcessing pre_pro(archive_filename);
PostProcessing post_pro(parse_data,pexport);
TLSGeom tls_geom(data, parse_data, pre_pro, post_pro,MPI_COMM_WORLD);
TLSSolver tls_solver(tls_geom, parse_data, pre_pro, post_pro);
QSFormulation formulation(tls_geom, tls_solver, data, parse_data, pre_pro, post_pro);
Algorithm algorithm(formulation, parse_data);
algorithm.run(pre_pro.getStep());
// wait every destruction done
MPI_Barrier(MPI_COMM_WORLD);
return eXlibris_tools::xMPIEnv::finalize();
}
......@@ -6,6 +6,8 @@
int main(int argc, char* argv[]) {
eXlibris_tools::xMPIEnv::init(argc,argv);
RegisterMaterials();
std::string data_filename, archive_filename;
......@@ -15,14 +17,18 @@ int main(int argc, char* argv[]) {
xParseData parse_data;
ReadData(data_filename, data, parse_data);
xfem::xExportGmshAscii pexport;
PreProcessing pre_pro(archive_filename);
PostProcessing post_pro(parse_data);
TLSGeom tls_geom(data, parse_data, pre_pro, post_pro);
PostProcessing post_pro(parse_data,pexport);
TLSGeom tls_geom(data, parse_data, pre_pro, post_pro,MPI_COMM_WORLD);
TLSSolver tls_solver(tls_geom, parse_data, pre_pro, post_pro);
QSFormulation formulation(tls_geom, tls_solver, data, parse_data, pre_pro, post_pro);
Algorithm algorithm(formulation, parse_data);
algorithm.run(pre_pro.getStep());
return 0;
// wait every destruction done
MPI_Barrier(MPI_COMM_WORLD);
return eXlibris_tools::xMPIEnv::finalize();
}
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