Commit bca90db0 authored by Kevin Moreau's avatar Kevin Moreau

Major update, lot of changes



git-svn-id: https://svn.ec-nantes.fr/eXlibris/Applis/DamageBandDyn@1751 fbbead7c-fb4d-4173-aa67-51132c73c120
parent 3b6eaca8
# ------------------------------------------------------------------------
# CMakeLists.txt for DamageBandDyn
# more documentation in http://www.cmake.org/cmake/help/cmake2.6docs.html
# ------------------------------------------------------------------------
#------------------------------------------------------------------------
cmake_minimum_required(VERSION 2.6)
project(DamageBandDyn)
......@@ -15,48 +15,49 @@ list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../../../Util/cmakeUti
include(common_functions)
include(devel_common_functions)
find_devroot(DEVROOT)
define_archos_suffixe(ARCHOS)
find_devroot(DEVROOT)
enable_testing()
# ---------------------------------------
# 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_XCRACK "use Xcrack" OFF)
option(USE_XTLS "use xTLS" ON)
option(USE_XCUT "use Xcut" ON)
option(USE_GEOM "use Geom" ON)
option(GZIPOUTPUT "use of boost iostream to compress output" OFF)
option(USE_OCTREE "use Octree" ON)
option(USE_MATLIB "use MatLib" ON)
if(USE_OCTREE OR USE_MATLIB)
set(USE_XINTERFACE ON)
endif(USE_OCTREE OR USE_MATLIB)
option(TIME_MONITORING "use time monitoring" ON)
option(USE_ANN "use library ANN" ON)
option(USE_CGAL "use library CGAL" OFF)
option(USE_CGAL "use library CGAL" ON)
if(USE_ANN OR USE_CGAL)
set(USE_GEOM ON)
endif(USE_ANN OR USE_CGAL)
option(USE_TR1 "use TR1" ON)
option(USE_MTL "use library MTL" ON)
option(USE_MTL "use library MTL" OFF)
# Generating extra stuffs
option(MAKE_DOC "make Doxygen" ON)
option(MAKE_TAGS "make Emacs tags" ON)
option(MAKE_TAGS "make tags" ON)
# Preprocessor variables
option(WITH_GZIPOUTPUT "use of boost iostream to compress output" ON)
option(TIME_MEASURE "use time monitoring" ON)
option(STATE_RESTORATION "use state restoration" ON)
option(EXTENDED_EXPORT "use extended export possibility" ON)
# --------------------
# treatment of options
# --------------------
set_internal_def()
set_external_def()
# ----------------------------------------
# Find external package not treated so far
# ----------------------------------------
find_and_set(Xfem XFEM DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
if(USE_TRELLIS)
add_definitions( -DHAVE_AOMD)
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)
......@@ -66,28 +67,40 @@ if(USE_XINTERFACE)
${SolverBase_INCLUDE_DIR}
${SOLVERINTERFACES_INCLUDE_DIR}
${SuperLu_INCLUDE_DIR}
${Taucs_INCLUDE_DIR}
)
endif(USE_XINTERFACE)
if (USE_XTLS)
if(USE_XCRACK)
add_definitions( -DHAVE_XCRACK)
find_and_set(Xcrack XCRACK DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_XCRACK)
if(USE_XTLS)
find_and_set(xTLS xTLS DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_XTLS)
if(USE_XCUT)
find_and_set(Xcut XCUT DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_XCUT)
if(USE_GEOM)
find_and_set(Geom GEOM DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_GEOM)
if(GZIPOUTPUT)
add_definitions( -DWITH_GIPZOUTPUT)
endif(GZIPOUTPUT)
if(USE_OCTREE)
add_definitions( -DHAVE_OCTREE)
find_and_set(Octree OCTREE DAMAGEBANDDYN_EXTERNAL_INCLUDES NOTHING)
endif(USE_OCTREE)
if(TIME_MONITORING)
add_definitions( -DTIME_MONITORING)# TODO in code pass from TIMING to TIME_MONITORING
endif(TIME_MONITORING)
if(USE_MATLIB)
add_definitions( -DHAVE_MATLIB)
find_package(MATLIB)
list(APPEND DAMAGEBANDDYN_EXTERNAL_INCLUDES ${MATLIB_INCLUDE_DIR} ${MatLibInterface_INCLUDE_DIR})
endif(USE_MATLIB)
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)
......@@ -102,6 +115,7 @@ if(USE_TR1)
endif(USE_TR1)
if(USE_MTL)
add_definitions( -DHAVE_MTL)
find_package(MTL)
list(APPEND DAMAGEBANDDYN_EXTERNAL_INCLUDES ${MTL_INCLUDE_DIR}/..)
endif(USE_MTL)
......@@ -113,23 +127,34 @@ if(MAKE_DOC)
COMMAND doxygen ${DEVROOT}/Applis/DamageBandDyn/doc/doxygen.config)
endif(MAKE_DOC)
if(WITH_GZIPOUTPUT)
add_definitions( -DWITH_GIPZOUTPUT)
endif(WITH_GZIPOUTPUT)
if(TIME_MEASURE)
add_definitions( -DTIME_MONITORING)
endif(TIME_MEASURE)
if(STATE_RESTORATION)
add_definitions( -DSTATE_RESTORATION)
endif(STATE_RESTORATION)
if(EXTENDED_EXPORT)
add_definitions( -DEXTENDED_EXPORT)
endif(EXTENDED_EXPORT)
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/*)
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)
add_library(DamageBandDyn ${BUILD_SHARED_LIBS} ${src_files})
include_directories(${DAMAGEBANDDYN_EXTERNAL_INCLUDES})
......@@ -138,8 +163,3 @@ set_target_properties(DamageBandDyn PROPERTIES COMPILE_FLAGS "-w -Wno-deprecated
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
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
// DamageBandDyn
#include "AsymmetricBehaviour.h"
#include "refactoring.h"
#ifdef PLANESTRESS
EvalPlaneStressAsymetricStress::EvalPlaneStressAsymetricStress(const xfem::xEval<double>& eval_strain_trace_, const xfem::xEval<xTensor2>& eval_strain_pos_part_, const xfem::xEval<xTensor2>& eval_strain_neg_part_, const xfem::xEval<double>& eval_damage_, const xfem::xEval<double>& eval_lambda_, const xfem::xEval<double>& eval_mu_) :
eval_strain_trace(eval_strain_trace_),
eval_strain_pos_part(eval_strain_pos_part_),
eval_strain_neg_part(eval_strain_neg_part_),
eval_damage(eval_damage_),
eval_lambda(eval_lambda_),
eval_mu(eval_mu_),
id(xfem::xVector(1.,0.,0.), xfem::xVector(0.,1.,0.), xfem::xVector(0.,0.,1.))
{
xfem::xTensorProduct tens_prod;
e3e3 = tens_prod(xfem::xVector(0.,0.,1.), xfem::xVector(0.,0.,1.));
}
void EvalPlaneStressAsymetricStress::operator()(const xfem::xGeomElem* geo_appro, const xfem::xGeomElem* geo_integ, result_type& res) const
{
double d, lambda, mu, tr, pos_part_tr, neg_part_tr;//, eps33;
xfem::xTensor2 pos_part, neg_part;
eval_strain_trace(geo_appro, geo_integ, tr);
eval_strain_pos_part(geo_appro, geo_integ, pos_part);
eval_strain_neg_part(geo_appro, geo_integ, neg_part);
eval_damage(geo_appro, geo_integ, d);
eval_lambda(geo_appro, geo_integ, lambda);
eval_mu(geo_appro, geo_integ, mu);
pos_part_tr = pospart(tr);
neg_part_tr = negpart(tr);
//eps33 = - (1.-d)*pos_part_tr - neg_part_tr;
res = pos_part*(2.*mu*(1.-d)) + neg_part*(2.*mu) + id*(lambda*( (1.-d)*pos_part_tr + neg_part_tr /*+ eps33*/ ) );
res(2,2) = 0.;
}
EvalPlaneStressAsymmetricStrain::EvalPlaneStressAsymmetricStrain(const xfem::xEval<xfem::xTensor2>& eval_strain_, const xfem::xEval<double>& eval_damage_) :
eval_strain(eval_strain_),
eval_damage(eval_damage_)
{}
void EvalPlaneStressAsymmetricStrain::operator()(const xfem::xGeomElem* geo_appro, const xfem::xGeomElem* geo_integ, result_type& res) const
{
double d, tr;
eval_strain(geo_appro, geo_integ, res);
eval_damage(geo_appro, geo_integ, d);
tr = res(0,0) + res(1,1);
res(2,2) = - (1.-d)*pospart(tr) - negpart(tr);
}
#endif
This diff is collapsed.
#ifndef AsymmetricBehaviour_imp_h_
#define AsymmetricBehaviour_imp_h_
namespace asymmetric_behaviour
{
template <typename TENSOR2>
EvalStress<TENSOR2>::
EvalStress(const xfem::xEval<result_type>& eval_strain_pos_part_,
const xfem::xEval<result_type>& eval_strain_neg_part_,
const xfem::xEval<double>& eval_damage_,
const xfem::xEval<double>& eval_lambda_,
const xfem::xEval<double>& eval_mu_,
DamageFct pos_part_fct_,
DamageFct neg_part_fct_) :
eval_strain_pos_part(eval_strain_pos_part_),
eval_strain_neg_part(eval_strain_neg_part_),
eval_damage(eval_damage_),
eval_lambda(eval_lambda_),
eval_mu(eval_mu_),
pos_part_fct(pos_part_fct_),
neg_part_fct(neg_part_fct_)
{}
template <typename TENSOR2>
void EvalStress<TENSOR2>::
operator()(const xfem::xGeomElem* geo_appro, const xfem::xGeomElem* geo_integ, result_type& res) const
{
eval_damage(geo_appro, geo_integ, buf);
#ifdef DEADZONESPECIALLAW
if (buf==1.) {
res=xfem::xTensor2(0.);
return;
}
#endif
eval_strain_pos_part(geo_appro, geo_integ, res);
eval_strain_neg_part(geo_appro, geo_integ, tensor);
if (buf==0.) {
res+=tensor;
tr=res.trace();
eval_mu(geo_appro, geo_integ, buf);
buf*=2.;
res*=buf;
eval_lambda(geo_appro, geo_integ, buf);
tr*=buf;
res(0,0)+=tr;
res(1,1)+=tr;
res(2,2)+=tr;
return;
}
tr=res.trace();
tr+=tensor.trace();
buf2=pos_part_fct(buf);
res*=buf2;
if (tr>0.) {
tr*=buf2;
}
buf2=neg_part_fct(buf);
tensor*=buf2;
if (tr<0.) {
tr*=buf2;
}
res+=tensor;
eval_mu(geo_appro, geo_integ, buf);
buf*=2.;
res*=buf;
if (tr==0.) {
return;
}
eval_lambda(geo_appro, geo_integ, buf);
tr*=buf;
res(0,0)+=tr;
res(1,1)+=tr;
res(2,2)+=tr;
}
template <typename TENSOR2>
EvalEnergyReleaseRate<TENSOR2>::
EvalEnergyReleaseRate(const xfem::xEval<Tensor2>& eval_strain_pos_part_,
const xfem::xEval<Tensor2>& eval_strain_neg_part_,
const xfem::xEval<double>& eval_damage_,
const xfem::xEval<double>& eval_lambda_,
const xfem::xEval<double>& eval_mu_,
DamageFctDerivative pos_part_fct_der_,
DamageFctDerivative neg_part_fct_der_) :
eval_strain_pos_part(eval_strain_pos_part_),
eval_strain_neg_part(eval_strain_neg_part_),
eval_damage(eval_damage_),
eval_lambda(eval_lambda_),
eval_mu(eval_mu_),
pos_part_fct_der(pos_part_fct_der_),
neg_part_fct_der(neg_part_fct_der_)
{}
template <typename TENSOR2>
void EvalEnergyReleaseRate<TENSOR2>::
operator()(const xfem::xGeomElem* geo_appro, const xfem::xGeomElem* geo_integ, result_type& res) const
{
eval_damage(geo_appro, geo_integ, buf);
buf2=-pos_part_fct_der(buf);
eval_strain_pos_part(geo_appro, geo_integ, tensor);
tr=tensor.trace();
res=tensor.contract(tensor);
res*=buf2;
buf=-neg_part_fct_der(buf);
eval_strain_neg_part(geo_appro, geo_integ, tensor);
res+=buf*tensor.contract(tensor);
tr+=tensor.trace();
if (tr<0.) {
tr*=tr;
tr*=buf;
}
else {
tr*=tr;
tr*=buf2;
}
eval_mu(geo_appro, geo_integ, buf);
res*=buf;
if (tr==0.) {
return;
}
eval_lambda(geo_appro, geo_integ, buf);
buf*=0.5;
tr*=buf;
res+=tr;
}
template <typename TENSOR2>
EvalStressRateBefore<TENSOR2>::
EvalStressRateBefore(const xfem::xEval<result_type>& eval_strain_pos_part_,
const xfem::xEval<result_type>& eval_strain_neg_part_,
const xfem::xEval<double>& eval_damage_,
const xfem::xEval<double>& eval_damage_rate_,
const xfem::xEval<double>& eval_lambda_,
const xfem::xEval<double>& eval_mu_,
DamageFctDerivative pos_part_fct_der_,
DamageFctDerivative neg_part_fct_der_) :
eval_strain_pos_part(eval_strain_pos_part_),
eval_strain_neg_part(eval_strain_neg_part_),
eval_damage(eval_damage_),
eval_damage_rate(eval_damage_rate_),
eval_lambda(eval_lambda_),
eval_mu(eval_mu_),
pos_part_fct_der(pos_part_fct_der_),
neg_part_fct_der(neg_part_fct_der_)
{}
template <typename TENSOR2>
void EvalStressRateBefore<TENSOR2>::
operator()(const xfem::xGeomElem* geo_appro, const xfem::xGeomElem* geo_integ, result_type& res) const
{
eval_damage_rate(geo_appro, geo_integ, d_dot);
if (d_dot==0.) {
res=xfem::xTensor2(0.);
return;
}
eval_damage(geo_appro, geo_integ, buf);
buf2 = pos_part_fct_der(buf);
buf = neg_part_fct_der(buf);
eval_strain_pos_part(geo_appro, geo_integ, res);
tr=res.trace();
res*=buf2;
eval_strain_neg_part(geo_appro, geo_integ, tensor);
tr+=tensor.trace();
tensor*=buf;
res+=tensor;
if (tr<0.) {
tr*=buf;
}
else {
tr*=buf2;
}
eval_mu(geo_appro, geo_integ, buf);
buf*=2.;
res*=(d_dot*buf);
if (tr==0.) {
return;
}
eval_lambda(geo_appro, geo_integ, buf);
tr*=(d_dot*buf);
res(0,0)+=tr;
res(1,1)+=tr;
res(2,2)+=tr;
}
template <typename TENSOR2>
EvalStressRateAfter<TENSOR2>::
EvalStressRateAfter(const xfem::xEval<result_type>& eval_strain_rate_pos_part_,
const xfem::xEval<result_type>& eval_strain_rate_neg_part_,
const xfem::xEval<double>& eval_damage_,
const xfem::xEval<double>& eval_lambda_,
const xfem::xEval<double>& eval_mu_,
DamageFct pos_part_fct_,
DamageFct neg_part_fct_) :
eval_strain_rate_pos_part(eval_strain_rate_pos_part_),
eval_strain_rate_neg_part(eval_strain_rate_neg_part_),
eval_damage(eval_damage_),
eval_lambda(eval_lambda_),
eval_mu(eval_mu_),
pos_part_fct(pos_part_fct_),
neg_part_fct(neg_part_fct_)
{}
template <typename TENSOR2>
void EvalStressRateAfter<TENSOR2>::
operator()(const xfem::xGeomElem* geo_appro, const xfem::xGeomElem* geo_integ, result_type& res) const
{
eval_damage(geo_appro, geo_integ, buf);
eval_strain_rate_pos_part(geo_appro, geo_integ, res);
eval_strain_rate_neg_part(geo_appro, geo_integ, tensor);
if (buf==0.) {
res+=tensor;
tr=res.trace();
eval_mu(geo_appro, geo_integ, buf);
buf*=2.;
res*=buf;
eval_lambda(geo_appro, geo_integ, buf);
tr*=buf;
res(0,0)+=tr;
res(1,1)+=tr;
res(2,2)+=tr;
return;
}
tr=res.trace();
tr+=tensor.trace();
buf2=pos_part_fct(buf);
res*=buf2;
if (tr>0.) {
tr*=buf2;
}
buf2=neg_part_fct(buf);
tensor*=buf2;
if (tr<0.) {
tr*=buf2;
}
res+=tensor;
eval_mu(geo_appro, geo_integ, buf);
buf*=2.;
res*=buf;
if (tr==0.) {
return;
}
eval_lambda(geo_appro, geo_integ, buf);
tr*=buf;
res(0,0)+=tr;
res(1,1)+=tr;
res(2,2)+=tr;
}
}
#endif
#include "CFL.h"
double ComputeMinInsideRadius::operator()(double val, AOMD::mEntity* e) const
{
Trellis_Util::mPoint p1 = static_cast<AOMD::mVertex*>(e->get(0,0))->point();
Trellis_Util::mPoint p2 = static_cast<AOMD::mVertex*>(e->get(0,1))->point();
Trellis_Util::mPoint p3 = static_cast<AOMD::mVertex*>(e->get(0,2))->point();
double a = Distance(p1, p2);
double b = Distance(p2, p3);
double c = Distance(p3, p1);
double k = 0.5*(a+b+c);
double h = std::sqrt( k*(k-a)*(k-b)*(k-c) ) / k;
return std::min(h, val);
}
double ComputeMaxInsideRadius::operator()(double val, AOMD::mEntity* e) const
{
Trellis_Util::mPoint p1 = static_cast<AOMD::mVertex*>(e->get(0,0))->point();
Trellis_Util::mPoint p2 = static_cast<AOMD::mVertex*>(e->get(0,1))->point();
Trellis_Util::mPoint p3 = static_cast<AOMD::mVertex*>(e->get(0,2))->point();
double a = Distance(p1, p2);
double b = Distance(p2, p3);
double c = Distance(p3, p1);
double k = 0.5*(a+b+c);
double h = std::sqrt( k*(k-a)*(k-b)*(k-c) ) / k;
return std::max(h, val);
}
std::pair<double,double> ComputeMinMaxInsideRadius::operator()(std::pair<double,double> val, AOMD::mEntity* e) const
{
Trellis_Util::mPoint p1 = static_cast<AOMD::mVertex*>(e->get(0,0))->point();
Trellis_Util::mPoint p2 = static_cast<AOMD::mVertex*>(e->get(0,1))->point();
Trellis_Util::mPoint p3 = static_cast<AOMD::mVertex*>(e->get(0,2))->point();
double a = Distance(p1, p2);
double b = Distance(p2, p3);
double c = Distance(p3, p1);
double k = 0.5*(a+b+c);
double h = std::sqrt( k*(k-a)*(k-b)*(k-c) ) / k;
return std::make_pair(std::min(h, val.first), std::max(h, val.second));
}
double ComputePamCrashLength::operator()(double val, AOMD::mEntity* e) const
{
double h_min = val;
xfem::xGeomElem geom_appro(e);
double meas_appro = geom_appro.getMeasure();
const int l = e->getLevel() - 1;
for (int i=0; i<e->size(l); ++i)
{
AOMD::mEntity* face = e->get(l,i);
xfem::xGeomElem geom_face(face);
double meas_face = geom_face.getMeasure();
double h = meas_appro/meas_face;
h_min = std::min(h, h_min);
}
return h_min;
}
ComputePamCrashTimeStep::ComputePamCrashTimeStep(const xfem::xEval<double>& eval_celerity_) :
eval_celerity(eval_celerity_)
{}
double ComputePamCrashTimeStep::operator()(double val, AOMD::mEntity* e) const
{
double dt_min = val;
xfem::xGeomElem geom_appro(e);
double celerity;
eval_celerity(&geom_appro, &geom_appro, celerity);
double meas_appro = geom_appro.getMeasure();
const int l = e->getLevel() - 1;
for (int i=0; i<e->size(l); ++i)
{
AOMD::mEntity* face = e->get(l,i);
xfem::xGeomElem geom_face(face);
double meas_face = geom_face.getMeasure();
double h = meas_appro/meas_face;
dt_min = std::min(h/celerity, dt_min);
}
return dt_min;
}
......@@ -5,123 +5,157 @@
*/
#ifndef _CFL_h
#define _CFL_h
#include <cmath>
// SolverInterface
#include "xLapackInterface.h"
// Xfem
#include "xAlgorithm.h"
// AOMD
#include "mPoint.h"
#include "mEntity.h"
#include "mVertex.h"
// std
#include <cmath>
// Xfem
#include "xGeomElem.h"
namespace AOMD
{
class mVertex;
}
namespace xfem
{
template <typename T>
class xEval;
}
inline double Distance(const Trellis_Util::mPoint& p1, const Trellis_Util::mPoint& p2)
/// Get distance between two points
inline double Distance(const Trellis_Util::mPoint& p1, const Trellis_Util::mPoint& p2)
{
return std::sqrt((p1(0)-p2(0))*(p1(0)-p2(0)) +
(p1(1)-p2(1))*(p1(1)-p2(1)) +
(p1(2)-p2(2))*(p1(2)-p2(2)));
return std::sqrt((p1(0)-p2(0))*(p1(0)-p2(0)) +
(p1(1)-p2(1))*(p1(1)-p2(1)) +
(p1(2)-p2(2))*(p1(2)-p2(2)));
}
/// Compute minimal inside radius of a triangle element, no check
struct ComputeMinInsideRadius
{
double operator()(double val, AOMD::mEntity* e)
{
Trellis_Util::mPoint p1 = ((AOMD::mVertex*) e->get(0,0))->point();
Trellis_Util::mPoint p2 = ((AOMD::mVertex*) e->get(0,1))->point();
Trellis_Util::mPoint p3 = ((AOMD::mVertex*) e->get(0,2))->point();
double a = Distance(p1, p2);
double b = Distance(p2, p3);
double c = Distance(p3, p1);
double k = 0.5*(a+b+c);
double h = std::sqrt( k*(k-a)*(k-b)*(k-c) ) / k;
return std::min(h, val);
}
double operator()(double val, AOMD::mEntity* e) const;
};
/// Compute maximal inside radius of a triangle element, no check
struct ComputeMaxInsideRadius
{
double operator()(double val, AOMD::mEntity* e)
{
Trellis_Util::mPoint p1 = ((AOMD::mVertex*) e->get(0,0))->point();
Trellis_Util::mPoint p2 = ((AOMD::mVertex*) e->get(0,1))->point();
Trellis_Util::mPoint p3 = ((AOMD::mVertex*) e->get(0,2))->point();
double a = Distance(