Commit a1e451a2 authored by Kevin Moreau's avatar Kevin Moreau

commit everything since merging of DamageBand and DamageBandDyn
through xTLS



git-svn-id: https://svn.ec-nantes.fr/eXlibris/Applis/DamageBandDyn@1556 fbbead7c-fb4d-4173-aa67-51132c73c120
parent 73649d31
# ------------------------------------------------------------------------
# CMakeLists.txt for DamageBandDyn
# more documentation in http://www.cmake.org/cmake/help/cmake2.6docs.html
# ------------------------------------------------------------------------
cmake_minimum_required(VERSION 2.6)
project(DamageBandDyn)
message(STATUS "-------------------------")
message(STATUS "PROJECT: DamageBandDyn --")
message(STATUS "-------------------------")
list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../../../Util/cmakeUtil/)
include(common_functions)
include(devel_common_functions)
define_archos_suffixe(ARCHOS)
find_devroot(DEVROOT)
# ---------------------------------------
# setting options
# option can be ON or OFF
# the test if() is true if it's ON
# ---------------------------------------
option(USE_TRELLIS "use Trellis" ON)
option(USE_OCTREE "use Octree" ON)
option(USE_XINTERFACE "use Xinterface" ON)
option(USE_XTLS "use xTLS" ON)
option(USE_GEOM "use Geom" ON)
option(GZIPOUTPUT "use of boost iostream to compress output" OFF)
option(TIME_MONITORING "use time monitoring" ON)
option(USE_ANN "use library ANN" ON)
option(USE_CGAL "use library CGAL" OFF)
option(USE_TR1 "use TR1" ON)
option(USE_MTL "use library MTL" ON)
option(MAKE_DOC "make Doxygen" ON)
option(MAKE_TAGS "make Emacs tags" ON)
# --------------------
# treatment of options
# --------------------
set_internal_def()
set_external_def()
# ----------------------------------------
# Find external package not treated so far
# ----------------------------------------
if(USE_TRELLIS)
find_and_set(Trellis TRELLIS DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_TRELLIS)
if(USE_OCTREE)
find_and_set(Octree OCTREE DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_OCTREE)
find_and_set(Xfem XFEM DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
if(USE_XINTERFACE)
find_package(Xinterfaces)
set_from_prefix(Xinterfaces DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
# temporary
list(APPEND DAMAGEBANDDYN_EXTERNAL_INCLUDES
${XINTERFACES_INCLUDE_DIR}
${SolverBase_INCLUDE_DIR}
${SOLVERINTERFACES_INCLUDE_DIR}
${SuperLu_INCLUDE_DIR}
)
endif(USE_XINTERFACE)
if (USE_XTLS)
find_and_set(xTLS xTLS DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_XTLS)
if(USE_GEOM)
find_and_set(Geom GEOM DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_GEOM)
if(GZIPOUTPUT)
add_definitions( -DWITH_GIPZOUTPUT)
endif(GZIPOUTPUT)
if(TIME_MONITORING)
add_definitions( -DTIME_MONITORING)# TODO in code pass from TIMING to TIME_MONITORING
endif(TIME_MONITORING)
if(USE_ANN)
add_definitions( -DHAVE_ANN)
#find_and_set(ANN ANN DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
find_package(ANN)
list(APPEND DAMAGEBANDDYN_EXTERNAL_INCLUDES ${ANN_INCLUDE_DIR})
endif(USE_ANN)
if(USE_CGAL)
add_definitions( -DHAVE_CGAL)
find_and_set(CGAL CGAL DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_CGAL)
if(USE_TR1)
find_package(TR1)
endif(USE_TR1)
if(USE_MTL)
find_package(MTL)
list(APPEND DAMAGEBANDDYN_EXTERNAL_INCLUDES ${MTL_INCLUDE_DIR}/..)
endif(USE_MTL)
if(MAKE_DOC)
add_custom_target(doc
COMMAND echo "Generating DamageBandDyn Doxygen"
WORKING_DIRECTORY ${DEVROOT}/Applis/DamageBandDyn/doc
COMMAND doxygen ${DEVROOT}/Applis/DamageBandDyn/doc/doxygen.config)
endif(MAKE_DOC)
set(BUILD_SHARED_LIBS "SHARED")
#set(LIBRARY_OUTPUT_PATH ${CMAKE_CURRENT_SOURCE_DIR}/lib/${ARCHOS})
#set(CMAKE_CURRENT_LIST_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src)
#find_package(BLAS2)
#find_package(LAPACK2)
#find_package(METIS)
#find_package(SolverInterfaces)
#find_package(Xinterfaces)
#find_package(SUPERLU)
#find_package(TAUCS)
#find_package(Boost)
add_definitions(-Wall)
file(GLOB src_files src/*)
add_library(DamageBandDyn ${BUILD_SHARED_LIBS} ${src_files})
include_directories(${DAMAGEBANDDYN_EXTERNAL_INCLUDES})
set_target_properties(DamageBandDyn PROPERTIES COMPILE_FLAGS "-w -Wno-deprecated ${CXXFLAGS}")
install(TARGETS DamageBandDyn DESTINATION ${DEVROOT}/Applis/DamageBandDyn/DamageBandDyn/lib/${ARCHOS})
if(MAKE_TAGS)
add_custom_target(tags
COMMAND echo "Generating DamageBandDyn Tags"
COMMAND etags --members --declarations ${src_files} --append --output ${DEVROOT}/Applis/DamageBandDyn/tags/TAGS)
endif(MAKE_TAGS)
\ No newline at end of file
#include "xAlgorithm.h"
#include "xData.h"
#include "Algorithm.h"
#include "Eval.h"
void TreatmentOfNatEnv(const xfem::xField& displacement_field, xfem::xAssembler& assembler, const xfem::xIntegrationRule& integrator, const xfem::xData& data, const xfem::xBoundary& groups, const double& time, const xfem::xEntityFilter filter)
{
if (data.mesh->dim()>1)
{
for (xfem::xPhysicalEnv::const_iterator it = data.PhysicalEnv->begin(); it != data.PhysicalEnv->end(); ++it)
{
const xfem::xEnv& env = *it;
if (env.Phys == "TRACTION_X" || env.Phys == "TRACTION_Y" || env.Phys == "TRACTION_Z" )
{
assert(env.Type == FIX);
xfem::xVector val;
if (env.Phys == "TRACTION_X") val(0) = env.getValue();
if (env.Phys == "TRACTION_Y") val(1) = env.getValue();
if (env.Phys == "TRACTION_Z") val(2) = env.getValue();
xfem::xEvalConstant<xfem::xVector> constant(val);
EvalSpaceFunctionCrossTimeFunction<xfem::xVector> traction(constant, env.getEvolution());
traction.setTime(time);
xfem::xFormLinearWithLoad<xfem::xValOperator<xfem::xIdentity<xfem::xVector> >, EvalSpaceFunctionCrossTimeFunction<xfem::xVector> > lin(traction);
xfem::xClassRegion bc(data.mesh, env.Entity, env.getDimension());
xfem::xFilteredRegion<xfem::xClassIter, xfem::xEntityFilter> fr_bc(bc.begin(), bc.end(), filter);
xfem::Assemble(lin, assembler, integrator, displacement_field, fr_bc.begin(), fr_bc.end(), xUpperAdjacency());
}
if (env.Phys == "PRESSURE")
{
assert(env.Type == FIX);
double val = env.getValue();
xEvalNormalbis eval_normal;
xfem::xEvalConstant<double> constant(val);
xfem::xEvalBinary<xfem::xMult<xfem::xVector, double, xfem::xVector> > pressure(eval_normal, constant);
EvalSpaceFunctionCrossTimeFunction<xfem::xVector> traction(pressure, env.getEvolution());
traction.setTime(time);
xfem::xFormLinearWithLoad<xfem::xValOperator<xfem::xIdentity<xfem::xVector> >, EvalSpaceFunctionCrossTimeFunction<xfem::xVector> > lin(traction);
xfem::xClassRegion bc(data.mesh, env.Entity, env.getDimension());
xfem::xFilteredRegion<xfem::xClassIter, xfem::xEntityFilter> fr_bc(bc.begin(), bc.end(), filter);
xfem::Assemble(lin, assembler, integrator, displacement_field, fr_bc.begin(), fr_bc.end(), xUpperAdjacency());
}
}
}
else
{
for (xfem::xPhysicalEnv::const_iterator it = data.PhysicalEnv->begin(); it != data.PhysicalEnv->end(); ++it)
{
const xfem::xEnv& env = *it;
if (env.Phys == "TRACTION_X")
{
assert(env.Type == FIX);
double val = 0.;
val = env.getValue();
xfem::xEvalConstant<double> constant(val);
EvalSpaceFunctionCrossTimeFunction<double> traction(constant, env.getEvolution());
traction.setTime(time);
xfem::xFormLinearWithLoad<xfem::xValOperator<xfem::xIdentity<double> >, EvalSpaceFunctionCrossTimeFunction<double> > lin(traction);
xfem::xClassRegion bc(data.mesh, env.Entity, env.getDimension());
xfem::xFilteredRegion<xClassIter, xEntityFilter> fr_bc(bc.begin(), bc.end(), filter);
xfem::Assemble(lin, assembler, integrator, displacement_field, fr_bc.begin(), fr_bc.end(), xUpperAdjacency());
}
}
}
}
void TreatmentOfEssEnv(const xfem::xField& displacement_field, const xfem::xField& velocity_field, const xfem::xField& acceleration_field, const xfem::xData& data, const double& time)
{
const int dim = data.mesh->dim();
for (xfem::xPhysicalEnv::const_iterator it = data.PhysicalEnv->begin(); it != data.PhysicalEnv->end(); ++it)
{
const xfem::xEnv& env = *it;
int geom = env.Geom;
if (geom==(dim-1))
{
std::string phys = env.Phys;
int type = env.Type;
int entity = env.Entity;
if (phys == "DISPLACEMENT_X" || phys == "DISPLACEMENT_Y" || phys == "DISPLACEMENT_Z")
{
assert(type==FIX);
double val = env.getValue();
// const xPieceWiseLinearDoubleToDouble& time_function = env.getEvolution();
// val *= time_function(time);
xfem::xClassRegion bc(data.mesh, entity, env.getDimension());
std::string kin("KINEMATIC_");
kin+=phys[phys.size()-1];
xfem::DirichletBoundaryCondition(displacement_field, kin, bc.begin(), bc.end(), val);
//DynamixDirichletBoundaryCondition(displacement_field, kin, bc.begin(), bc.end(), mass, val);
}
if (phys == "VELOCITY_X" || phys == "VELOCITY_Y" || phys == "VELOCITY_Z")
{
assert(type==FIX);
double val = env.getValue();
// const xPieceWiseLinearDoubleToDouble& time_function = env.getEvolution();
// val *= time_function(time);
xfem::xClassRegion bc(data.mesh, entity, env.getDimension());
std::string kin("KINEMATIC_");
kin+=phys[phys.size()-1];
xfem::DirichletBoundaryCondition(velocity_field, kin, bc.begin(), bc.end(), val);
}
if (phys == "ACCELERATION_X" || phys == "ACCELERATION_Y" || phys == "ACCELERATION_Z" )
{
assert(type==FIX);
double val = env.getValue();
// const xPieceWiseLinearDoubleToDouble& time_function = env.getEvolution();
// val *= time_function(time);
xfem::xClassRegion bc(data.mesh, entity, env.getDimension());
std::string kin("KINEMATIC_");
kin+=phys[phys.size()-1];
xfem::DirichletBoundaryCondition(acceleration_field, kin, bc.begin(), bc.end(), val);
}
}
}
}
#ifndef _Algorithm_h
#define _Algorithm_h
#include "xEntityFilter.h"
namespace xfem
{
class xAssembler;
class xData;
class xField;
class xIntegrationRule;
class xBoundary;
class xLevelSet;
}
/// \brief Hacked version of TreatmentOfNatEnv for dynamics (to have
/// 1D and spatial partitioning)
void TreatmentOfNatEnv(const xfem::xField& displacement_field, xfem::xAssembler& assembler, const xfem::xIntegrationRule& integrator, const xfem::xData& data, const xfem::xBoundary& groups, const double& time, const xfem::xEntityFilter filter=xfem::xAcceptAll());
/// \brief Dirichlet boundary conditions in a dynamics context
void TreatmentOfEssEnv(const xfem::xField& displacement_field, const xfem::xField& velocity_field, const xfem::xField& acceleration_field, const xfem::xData& data, const double& time);
/// \brief attach specific BC to a vertex mEntity*
template <typename EVAL, typename ITER>
void SpecificDirichletBoundaryCondition(const EVAL& eval_, const ITER& begin_, const ITER& end_, std::map<AOMD::mEntity*, double>& front_velocity_map_, const xfem::xEntityToEntity& entity_to_entity_);
/// \brief Initial conditions from L2Projection
template <typename ITER>
void ImposeInitialConditions(const xfem::xData& data_, xfem::xField& displacement_field_, xfem::xField& velocity_field_, const xfem::xIntegrationRule& integration_rule_, const ITER& begin_, const ITER& end_);
/// \brief Search for intitiation
template <typename EVAL1, typename EVAL2, typename ITER>
void SearchForInitiation(const EVAL1& eval_, const EVAL2& eval_threshold_, const xfem::xIntegrationRule& integration_rule_, const ITER& begin_, const ITER& end_, const double max_element_size_, xfem::xLevelSet& ls_in_, bool& first_initiated_, bool& newly_initiated_);
#include "Algorithm_imp.h"
#endif
#ifndef _Algorithm_imp_h
#define _Algorithm_imp_h
#ifdef DEBUG
#include <iostream>
#endif
// Xfem
#include "xAlgorithm.h"
#include "xAssembler.h"
#include "xData.h"
#include "xSimpleGeometry.h"
// SolverBase
#include "xLinearSystemSolverSuperLU.h"
// DamageBandDyn
#include "CommandOnGeomElem.h"
#include "EntityToEntity.h"
#include "Eval.h"
#include "LevelSetOperators.h"
template <class EVAL, class ITER>
void SpecificDirichletBoundaryCondition(const EVAL& eval_, const ITER& begin_, const ITER& end_, std::map<AOMD::mEntity*, double>& front_velocity_map_, const xfem::xEntityToEntity& entity_to_entity_)
{
for (ITER iter=begin_; iter!=end_; ++iter)
{
AOMD::mVertex* v = static_cast<AOMD::mVertex*>(*iter);
Trellis_Util::mPoint p = v->point();
AOMD::mEntity* entity = entity_to_entity_(v);
if (entity)
{
xfem::xGeomElem geom_elem(entity);
geom_elem.setUVWForXYZ(p);
double val;
eval_(&geom_elem, &geom_elem, val);
front_velocity_map_.insert(std::pair<AOMD::mEntity*, double>(v, val));
}
else
{
#ifdef DEBUG
std::cerr << "no entity to entity found in SpecificDirichletBoundaryCondition" << std::endl;
#endif
throw;
}
}
}
template <class ITER>
void ImposeInitialConditions(const xfem::xData& data_, xfem::xField& displacement_field_, xfem::xField& velocity_field_, const xfem::xIntegrationRule& integration_rule_, const ITER& begin_, const ITER& end_)
{
const int dim = data_.mesh->dim();
if (dim>1)
{
xfem::xVector initial_displacement(0.,0.,0.);
xfem::xVector initial_velocity(0.,0.,0.);
for (xfem::xPhysicalEnv::const_iterator it=data_.PhysicalEnv->begin(); it!=data_.PhysicalEnv->end(); ++it)
{
const xfem::xEnv& env = *it;
int geom = env.Geom;
if (geom==dim) //BC_SURFACE in 2D and BC_VOLUME in 3D
{
std::string phys = env.Phys;
double value = env.Val_fix;
if (phys == "DISPLACEMENT_X") initial_displacement(0) = value;
if (phys == "DISPLACEMENT_Y") initial_displacement(1) = value;
if (phys == "DISPLACEMENT_Z") initial_displacement(2) = value;
if (phys == "VELOCITY_X") initial_velocity(0) = value;
if (phys == "VELOCITY_Y") initial_velocity(1) = value;
if (phys == "VELOCITY_Z") initial_velocity(2) = value;
}
}
xfem::xEvalConstant<xVector> eval_initial_displacement(initial_displacement);
xfem::xEvalConstant<xVector> eval_initial_velocity(initial_velocity);
xfem::xAssemblerBasic<> assembler;
lalg::xLinearSystemSolverSuperLU<> solver_disp_lu;
L2Projection(displacement_field_, eval_initial_displacement, assembler, integration_rule_, solver_disp_lu, begin_, end_);
lalg::xLinearSystemSolverSuperLU<> solver_speed_lu;
L2Projection(velocity_field_, eval_initial_velocity, assembler, integration_rule_, solver_speed_lu, begin_, end_);
}
else
{
double initial_displacement = 0.;
double initial_velocity = 0.;
for (xfem::xPhysicalEnv::const_iterator it=data_.PhysicalEnv->begin(); it!=data_.PhysicalEnv->end(); ++it)
{
const xfem::xEnv& env = *it;
int geom = env.Geom;
if (geom==1) //BC_Line in 1D
{
std::string phys = env.Phys;
double value = env.Val_fix;
if (phys == "DISPLACEMENT_X") initial_displacement = value;
if (phys == "VELOCITY_X") initial_velocity = value;
}
}
xfem::xEvalConstant<double> eval_initial_displacement(initial_displacement);
xfem::xEvalConstant<double> eval_initial_velocity(initial_velocity);
xfem::xAssemblerBasic<> assembler;
lalg::xLinearSystemSolverSuperLU<> solver_disp_lu;
L2Projection(displacement_field_, eval_initial_displacement, assembler, integration_rule_, solver_disp_lu, begin_, end_);
lalg::xLinearSystemSolverSuperLU<> solver_speed_lu;
L2Projection(velocity_field_, eval_initial_velocity, assembler, integration_rule_, solver_speed_lu, begin_, end_);
}
}
template <typename EVAL1, typename EVAL2, typename ITER>
void SearchForInitiation(const EVAL1& eval_, const EVAL2& eval_threshold_, const xfem::xIntegrationRule& integration_rule_, const ITER& begin_, const ITER& end_, const double max_element_size_, xfem::xLevelSet& ls_in_, bool& first_initiated_, bool& newly_initiated_)
{
newly_initiated_ = false;
std::vector<Trellis_Util::mPoint> hot_points;
LocateHotPointsCommand<EVAL1, EVAL2> locate_hot_points(eval_, eval_threshold_, hot_points);
ApplyCommandOnIntegrationRule(locate_hot_points, integration_rule_, begin_, end_);
if (!hot_points.empty())
{
std::vector<xfem::xPointToDouble> ptd_vector;
for (std::vector<Trellis_Util::mPoint>::const_iterator it=hot_points.begin(); it!=hot_points.end(); ++it)
{
ptd_vector.push_back(xfem::xCompl(xfem::xSphere(*it, 4*max_element_size_)));
}
//ptd_vector.push_back(xfem::xCompl(xfem::xSphere(*hot_points.begin(), 4*max_element_size_)));
xfem::xInterVector inter_vector(ptd_vector);
//newly_initiated_ = true;
if (!first_initiated_) // initiate for the first time
{
// tmp
newly_initiated_ = true;
// tmp
first_initiated_ = true;
ls_in_.clear();
ls_in_.load(inter_vector);
message(__func__, "found first hot points");
}
else
{
//AddPointToDoubleToLevelSet add_point_to_double_to_level_set(inter_vector);
//add_point_to_double_to_level_set.visit(ls_in_, ls_in_.getSupport());
//message(__func__, "found new hot points");
}
}
}
#endif
#ifndef _CFL_h
#define _CFL_h
#include <cmath>
// AOMD
#include "mPoint.h"
#include "mEntity.h"