Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65464893
compute_arlequin_weight.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Jun 4, 00:02
Size
9 KB
Mime Type
text/x-c++
Expires
Thu, Jun 6, 00:02 (2 d)
Engine
blob
Format
Raw Data
Handle
18076299
Attached To
rLIBMULTISCALE LibMultiScale
compute_arlequin_weight.cc
View Options
/**
* @file compute_arlequin_weight.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "compute_arlequin_weight.hh"
#include "lib_continuum.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "ref_point_data.hh"
#include "trace_atom.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
///* --------------------------------------------------------------------------
///*/
// enum FunctionOrder { LINEAR = 1, THIRD = 3, CUSTOM = 4 };
///* -------------------------------------------------------------------------
///*/
//
// struct Weight1D {
// template <FunctionOrder order> static Real weight(Real x, Real R);
//};
/* ------------------------------------------------------------------------ */
Real
ComputeArlequinWeight
::
computeWeight
(
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>
// inline Weighting<FLAG, Dim, order>::Weighting(LMID geom_id) {
// this->geom = GeometryManager::getManager().getGeometry(geom_id);
//}
/* ------------------------------------------------------------------------ */
// template <UInt FLAG, UInt Dim, FunctionOrder order>
// template <typename T> inline Real ComputeArlequinWeight::computeWeight(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>
// inline Real Weighting<FLAG, Dim, order>::computeWeight(Real x, Real R) {
// return Weight1D::weight<order>(x, R);
//}
//
//
/* --------------------------------------------------------------------------
*/
template
<
UInt
Dim
>
inline
Real
ComputeArlequinWeight
::
weight
(
Vector
<
Dim
>
X
)
{
auto
&&
geom
=
GeometryManager
::
getGeometry
(
this
->
geom
);
switch
(
geom
->
getType
())
{
case
Geometry
::
BALL:
return
this
->
weight
(
X
,
*
std
::
dynamic_pointer_cast
<
Ball
>
(
geom
));
return
0.
;
case
Geometry
::
CUBE:
return
this
->
weight
(
X
,
*
std
::
dynamic_pointer_cast
<
Cube
>
(
geom
));
default
:
LM_TOIMPLEMENT
;
}
}
/* ------------------------------------------------------------------------ */
template
<
UInt
Dim
>
inline
Real
ComputeArlequinWeight
::
weight
(
Vector
<
Dim
>
X
,
Ball
&
g
)
{
Real
result
=
0
;
DUMP
(
"employed gemometry "
<<
g
,
DBG_DETAIL
);
auto
x
=
X
-
g
.
getCenter
().
cast
<
Real
>
().
block
<
Dim
,
1
>
(
0
,
0
);
Real
d
=
sqrt
(
x
.
dot
(
x
));
DUMP
(
"after correction x = "
<<
x
[
0
]
<<
" y = "
<<
x
[
1
]
<<
" z = "
<<
x
[
2
]
<<
" d = "
<<
d
,
DBG_DETAIL
);
result
=
computeWeight
(
d
-
g
.
rMin
(),
g
.
rMax
()
-
g
.
rMin
());
return
result
;
}
/* ------------------------------------------------------------------------ */
template
<
UInt
Dim
>
inline
Real
ComputeArlequinWeight
::
weight
(
Vector
<
Dim
>
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
;
}
/* ------------------------------------------------------------------------ */
ComputeArlequinWeight
::
ComputeArlequinWeight
(
const
LMID
&
id
)
:
LMObject
(
id
)
{
this
->
createInput
(
"input"
);
this
->
createOutput
(
"weight"
);
this
->
createArrayOutput
(
"weight"
);
}
/* --------------------------------------------------------------------------
*/
ComputeArlequinWeight
::~
ComputeArlequinWeight
()
{}
/* ------------------------------------------------------------------------ */
template
<
typename
Cont
>
enable_if_mesh
<
Cont
>
ComputeArlequinWeight
::
build
(
Cont
&
meshList
)
{
constexpr
UInt
Dim
=
Cont
::
Dim
;
auto
&
weightFE
=
this
->
getOutputAsArray
(
"weight"
);
UInt
nb
=
meshList
.
size
();
weightFE
.
assign
(
nb
,
0
);
if
(
!
nb
)
{
DUMP
(
"There are no nodes in the bridging"
,
DBG_WARNING
);
return
;
}
#ifndef LM_OPTIMIZED
auto
&
elemList
=
meshList
.
getContainerElems
();
DUMP
(
"We found "
<<
nb
<<
" nodes concerned in "
<<
elemList
.
size
()
<<
" elements"
,
DBG_INFO
);
#endif
// LM_OPTIMIZED
for
(
auto
&&
[
n
,
w
]
:
zip
(
meshList
,
weightFE
))
{
Vector
<
Dim
>
pos0
=
n
.
position0
();
w
=
this
->
weight
(
pos0
);
if
(
w
>
1.
)
w
=
1.
;
if
(
w
<
ZERO_LIMIT
)
w
=
quality
;
DUMP
(
"w = "
<<
w
,
DBG_ALL
);
}
}
/* ------------------------------------------------------------------------ */
template
<
typename
Cont
>
enable_if_point
<
Cont
>
ComputeArlequinWeight
::
build
(
Cont
&
pointList
)
{
constexpr
UInt
Dim
=
Cont
::
Dim
;
UInt
nb
=
pointList
.
size
();
auto
&
weightMD
=
this
->
getOutputAsArray
(
"weight"
);
weightMD
.
assign
(
nb
,
0
);
if
(
!
nb
)
{
DUMP
(
"We found no atoms in the bridging zone"
,
DBG_WARNING
);
return
;
}
#ifndef LM_OPTIMIZED
DUMP
(
"We found "
<<
nb
<<
" concerned atoms"
,
DBG_INFO
);
#endif
// LM_OPTIMIZED
for
(
auto
&&
[
at
,
w
]
:
zip
(
pointList
,
weightMD
))
{
Vector
<
Dim
>
pos0
=
at
.
position0
();
w
=
1.
-
this
->
weight
(
pos0
);
if
(
w
<
ZERO_LIMIT
)
w
=
quality
;
DUMP
(
"w = "
<<
w
,
DBG_ALL
);
}
}
/* --------------------------------------------------------------------------
*/
/* LMDESC ARLEQUINWEIGHT
This computes a typical weight function, usable for arlequin
coupling. Beware there are strong geometry assumptions (like syymetry).
*/
/* LMEXAMPLE
COMPUTE weight ARLEQUIN_WEIGHT INPUT md GEOM some_geometry
*/
inline
void
ComputeArlequinWeight
::
declareParams
()
{
ComputeInterface
::
declareParams
();
/* LMKEYWORD QUALITY
Because of the strong sense brought in the formulation
of the Lagrange constraints zero weights are prohibited.
Therefore the so-called quality factor defines the
replacement value for these zeros.
More details on the impact of the factor can be found in
*Ghost force reduction and spectral analysis of the 1D bridging
method*.
**Guillaume Anciaux, Olivier Coulaud, Jean Roman, Gilles Zerah.**
`(Link) <http://hal.inria.fr/inria-00300603/en/>`_
*/
this
->
parseKeyword
(
"QUALITY"
,
quality
);
/* LMKEYWORD GEOMETRY
Set the bridging/overlaping zone where the Lagrange multipliers are
to be computed.
*/
this
->
parseKeyword
(
"GEOMETRY"
,
geom
);
}
/* --------------------------------------------------------------------------
*/
DECLARE_COMPUTE_MAKE_CALL
(
ComputeArlequinWeight
)
__END_LIBMULTISCALE__
Event Timeline
Log In to Comment