Page MenuHomec4science
Files F82597777

weighting.hh
No OneTemporary
Actions

File Metadata

Created
Thu, Sep 12, 08:25

weighting.hh
View Options

/**
* @file weighting.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Nov 25 12:23:57 2013
*
* @brief This is a convenient class used to compute the weights associated with most methods based on a bridging zone
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_WEIGHTING_HH__
#define __LIBMULTISCALE_WEIGHTING_HH__
/* -------------------------------------------------------------------------- */
#define CONTINUFLAG 1
#define ATOMEFLAG 2
/* -------------------------------------------------------------------------- */
#include "geom_inter.hh"
#include "cube.hh"
#include "ball.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/* -------------------------------------------------------------------------- */
enum FunctionOrder{
LINEAR = 1,
THIRD = 3,
CUSTOM = 4
};
/* ------------------------------------------------------------------------- */
struct Weight1D {
template<FunctionOrder order>
static Real weight(Real x,Real R);
};
/* -------------------------------------------------------------------------- */
template <>
inline Real Weight1D::weight<LINEAR>(Real x,Real R){
return (x/R);
}
/* -------------------------------------------------------------------------- */
template <>
inline Real Weight1D::weight<THIRD>(Real x,Real R){
Real R3 = R*R*R;
return (-1.*x*x*(2.*x-3*R)/R3);
}
/* -------------------------------------------------------------------------- */
template <>
inline Real Weight1D::weight<CUSTOM>(Real x,Real R){
/* Real r=R,q=8*R/10,epsilon=0.1; */
/* Real a = -(epsilon*r+q-r)/(q-r)/(q-r)/q; */
/* Real b = (q*q*epsilon -r*r +q*q +epsilon*r*r)/(q-r)/(q-r)/q; */
/* Real c = - (q+epsilon*q-r)*r/(q-r)/(q-r); */
/* if(x < q) */
/* return (1-epsilon)/q*x; */
/* return (1-a*x*x - b*x -c); */
/* ********************** autre experience **************/
/* Real epsilon = .1; */
/* // Real q = R/16;//marche pas */
/* Real q = R/32; */
/* if (x < q) */
/* return 0; */
/* return 1/(R-q)*(x-q); */
return (x/R);
}
/* -------------------------------------------------------------------------- */
template <UInt FLAG, UInt Dim, FunctionOrder order>
class Weighting {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Weighting(GeomID geom_id);
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
//! compute a weight associated with coordinates
Real weight(Real * X);
protected:
//! specializations
Real weight(Real * X,Cube & c);
Real weight(Real * X,Ball & b);
//! compute a weight associated with 1D coordinates
Real computeWeight(Real x,Real R);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
Geometry * geom;
};
/* -------------------------------------------------------------------------- */
template <UInt FLAG, UInt Dim, FunctionOrder order>
inline Weighting<FLAG,Dim,order>::Weighting(GeomID geom_id){
this->geom = GeometryManager::getManager().getGeometry(geom_id);
}
/* -------------------------------------------------------------------------- */
template <UInt FLAG, UInt Dim, FunctionOrder order>
inline Real Weighting<FLAG,Dim,order>::weight(Real * X){
Real result=-1;
if (geom == NULL){
std::stringstream sstr;
sstr << "weighting of point at coordinate ";
for (UInt k = 0; k < Dim; ++k) sstr << X[k] << " ";
LM_FATAL(sstr.str()
<< " impossible since passed geometry is null");
}
Geometry & g = *geom;
switch(g.getType()){
case Geometry::CUBE : result = weight(X,dynamic_cast<Cube &>(*geom)) ; break;
case Geometry::BALL : result = weight(X,dynamic_cast<Ball &>(*geom)) ; break;
default: LM_TOIMPLEMENT;
}
if (result > 1.0)
result = 1.0;
if (result < 0.0)
result = 0.0;
if(FLAG==ATOMEFLAG) return (1.0 - result);
return result;
}
/* -------------------------------------------------------------------------- */
template <UInt FLAG, UInt Dim, FunctionOrder order>
inline Real Weighting<FLAG,Dim,order>::weight(Real * X,Cube & c){
Real result = 0;
if (c.getXmax()[1] < 0){
Real R = c.getXmax()[1] - c.getXmin()[1];
result = computeWeight(c.getXmax()[1]-X[1],R);
}
else if (c.getXmin()[1] > 0){
Real R = c.getXmax()[1] - c.getXmin()[1];
result = computeWeight(X[1]-c.getXmin()[1],R);
}
else if (c.getXmax()[2] < 0){
Real R = c.getXmax()[2] - c.getXmin()[2];
result = computeWeight(c.getXmax()[2]-X[2],R);
}
else if (c.getXmin()[2] > 0){
Real R = c.getXmax()[2] - c.getXmin()[2];
result = computeWeight(X[2]-c.getXmin()[2],R);
}
else if (c.getXmax()[0] < 0){
Real R = c.getXmax()[0] - c.getXmin()[0];
result = computeWeight(c.getXmax()[0]-X[0],R);
}
else if (c.getXmin()[0] > 0){
Real R = c.getXmax()[0] - c.getXmin()[0];
result = computeWeight(X[0]-c.getXmin()[0],R);
}
else
LM_FATAL("bridged zone can not be of size 0");
return result;
}
/* -------------------------------------------------------------------------- */
template <UInt FLAG, UInt Dim, FunctionOrder order>
inline Real Weighting<FLAG,Dim,order>::weight(Real * X, Ball & g){
Real x[3] = {0,0,0};
DUMP("la geometrie qui arrive ici " << g,DBG_DETAIL);
for (UInt k = 0; k < Dim; ++k)
x[k] = X[k] - g.getCenter(k);
Real d = sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
DUMP("after correction x = " << x[0] << " y = " << x[1] << " z = " << x[2] << " d = " << d,DBG_DETAIL);
return computeWeight(d-g.rMin(),g.rMax()-g.rMin());
}
/* -------------------------------------------------------------------------- */
template <UInt FLAG, UInt Dim, FunctionOrder order>
inline Real Weighting<FLAG,Dim,order>::computeWeight(Real x,Real R){
return Weight1D::weight<order>(x,R);
}
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_WEIGHTING_HH__ */

Event Timeline

c4science · Help