Homec4science
Diffusion LibMultiScale (features/reshape_compute_python) Browse

libmultiscale/src/coupling/weighting.hh
be6b3002c617features/reshape_comput…

Last Change
Raw File

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 "ball.hh"
#include "cube.hh"
#include "geom_inter.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
template <typename V> Real weight(V X);
protected:
//! specializations
template <typename V> Real weight(V X, Cube &c);
template <typename V> Real weight(V 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>
template <typename T>
inline Real Weighting<FLAG, Dim, order>::weight(T 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>
template <typename V>
inline Real Weighting<FLAG, Dim, order>::weight(V 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>
template <typename V>
inline Real Weighting<FLAG, Dim, order>::weight(V 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__ */
c4science · Help