Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F83662000
domain_paradis.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
Wed, Sep 18, 09:03
Size
32 KB
Mime Type
text/x-c
Expires
Fri, Sep 20, 09:03 (2 d)
Engine
blob
Format
Raw Data
Handle
20874959
Attached To
rLIBMULTISCALE LibMultiScale
domain_paradis.cc
View Options
/**
* @file domain_paradis.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Moseley Philip Arthur <philip.moseley@epfl.ch>
* @author Jaehyun Cho <jaehyun.cho@epfl.ch>
* @author Max Hodapp <max.hodapp@epfl.ch>
*
* @date Fri Jul 11 15:47:44 2014
*
* @brief ParaDiS domain
*
* @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/>.
*
*/
/* TODO:
* - the function addHybSegments() is currently a rudimentary implementation
* for simulations containing only one dislocations
* ==> method has to be generalized and additional checks (e.g. if hybrid
* segments are actually self-consistent, etc...) have to be performed to
* improve robustness
*
*/
/* -------------------------------------------------------------------------- */
#define TIMER
/* -------------------------------------------------------------------------- */
#include "lm_common.hh"
#define PARALLEL
#include "communicator.hh"
#include "domain_paradis.hh"
#include "geometry_manager.hh"
#include "stdio.h"
/* -------------------------------------------------------------------------- */
// ParaDiS includes
extern
"C"
{
// void TimerClearAll(Home_t *home);
#include "Home.h"
#include "Param.h"
#include "Typedefs.h"
/* -------------------------------------------------------------------------- */
Home_t
*
InitHome
();
void
ParadisStep
(
Home_t
*
home
);
void
RemeshAfterCoupling
(
Home_t
*
home
);
void
PointDisplacements
(
Home_t
*
home
);
void
SegmentDisplacementsOfAtom
(
libmultiscale
::
Vector
<
3u
>
&
A
,
libmultiscale
::
Vector
<
3u
>
&
B
,
libmultiscale
::
Vector
<
3u
>
&
C
,
libmultiscale
::
Vector
<
3u
>
&
burg
,
libmultiscale
::
Vector
<
3u
>
&
dispPointU
,
libmultiscale
::
Vector
<
3u
>
&
dispPointX
,
libmultiscale
::
Real
pois
);
void
ComputeClosurePoint
(
libmultiscale
::
Vector
<
3u
>
&
A
,
libmultiscale
::
Vector
<
3u
>
&
B
,
libmultiscale
::
Vector
<
3u
>
&
burg
,
libmultiscale
::
Vector
<
3u
>
&
C
);
void
ComputeProjectedPoint
(
libmultiscale
::
Vector
<
3u
>
&
A
,
libmultiscale
::
Vector
<
3u
>
&
C
,
libmultiscale
::
Vector
<
3u
>
&
burg
,
libmultiscale
::
Vector
<
3u
>
&
P
);
void
GetFieldPointStress
(
Home_t
*
home
,
real8
x
,
real8
y
,
real8
z
,
real8
totStress
[
3
][
3
]);
#ifdef FPES_ON
#include <fpcontrol.h>
#endif
}
#undef PARALLEL
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
void
DomainParaDiS
::
remesh
()
{
RemeshAfterCoupling
(
this
->
paradis_ptr
);
// //ParadisStep.c from line 336 to 418
// InitTopologyExemptions(this->paradis_ptr);
// ClearOpList(this->paradis_ptr);
// SortNodesForCollision(this->paradis_ptr);
// SplitMultiNodes(this->paradis_ptr);
// CrossSlip(this->paradis_ptr);
// HandleCollisions(this->paradis_ptr);
// ProximityCollisions(this->paradis_ptr);
// #ifdef PARALLEL
// #ifdef PARADIS_IN_LIBMULTISCALE
// CommSendRemesh(this->paradis_ptr);
// #endif
// #endif
// FixRemesh(this->paradis_ptr);
// Remesh(this->paradis_ptr);
// Migrate(this->paradis_ptr);
// #ifdef PARALLEL
// #ifdef PARADIS_IN_LIBMULTISCALE
// RecycleGhostNodes(this->paradis_ptr);
// #endif
// #endif
// SortNativeNodes(this->paradis_ptr);
// #ifdef PARALLEL
// #ifdef PARADIS_IN_LIBMULTISCALE
// CommSendGhosts(this->paradis_ptr);
// #endif
// #endif
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
computeDisplacements
(
std
::
vector
<
Real
>
&
positions
,
std
::
vector
<
Real
>
&
displacements
)
{
if
(
positions
.
empty
())
{
displacements
.
clear
();
return
;
}
if
(
!
this
->
paradis_ptr
)
{
return
;
}
// Pass data to Paradis.
unsigned
int
size
=
positions
.
size
();
LM_ASSERT
(
size
%
3
==
0
,
"paradis computeDisplacements requires 3d."
);
if
(
displacements
.
size
()
!=
size
)
displacements
.
resize
(
size
);
this
->
paradis_ptr
->
dispPointCount
=
size
/
3
;
this
->
paradis_ptr
->
dispPointX
=
&
positions
[
0
];
this
->
paradis_ptr
->
dispPointU
=
&
displacements
[
0
];
// Calculate the displacements.
PointDisplacements
(
this
->
paradis_ptr
);
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
computeBarnett
(
Vector
<
3u
>
&
A
,
Vector
<
3u
>
&
B
,
Vector
<
3u
>
&
C
,
Vector
<
3u
>
&
burg
,
Vector
<
3u
>
&
dispPointU
,
Vector
<
3u
>
&
dispPointX
,
Real
pois
)
{
SegmentDisplacementsOfAtom
(
A
,
B
,
C
,
burg
,
dispPointU
,
dispPointX
,
pois
);
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
computeClosurePoint
(
Vector
<
3u
>
&
A
,
Vector
<
3u
>
&
B
,
Vector
<
3u
>
&
burg
,
Vector
<
3u
>
&
C
)
{
ComputeClosurePoint
(
A
,
B
,
burg
,
C
);
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
computeProjectedPoint
(
Vector
<
3u
>
&
A
,
Vector
<
3u
>
&
C
,
Vector
<
3u
>
&
burg
,
Vector
<
3u
>
&
P
)
{
ComputeProjectedPoint
(
A
,
C
,
burg
,
P
);
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
init
()
{
auto
&
group
=
this
->
getCommGroup
();
if
(
!
group
.
amIinGroup
())
return
;
const
UInt
nbprocs
=
group
.
size
();
// Mostly taken from ParaDiS source files Main.c and ParadisInit.c
std
::
string
dlb_flag
=
"-r"
;
std
::
string
datafile_flag
=
"-d"
;
std
::
string
ctrlfile_flag
=
"-c"
;
std
::
string
fake_argv0
=
"IGNORED"
;
std
::
string
dlb_num
;
std
::
stringstream
(
dlb_num
)
<<
numDLBCycles
;
std
::
vector
<
const
char
*>
argv
({
fake_argv0
.
c_str
(),
datafile_flag
.
c_str
(),
data_filename
.
c_str
(),
dlb_flag
.
c_str
(),
dlb_num
.
c_str
(),
ctrlfile_flag
.
c_str
(),
control_filename
.
c_str
()});
if
(
mobility_filename
!=
""
)
argv
.
push_back
(
mobility_filename
.
c_str
());
unsigned
int
argc
=
argv
.
size
();
/*
* On some systems, the getrusage() call made by Meminfo() to get
* the memory resident set size does not work properly. In those
* cases, the function will try to return the current heap size
* instead. This initial call allows meminfo() to get a copy of
* the original heap pointer so subsequent calls can calculate the
* heap size by taking the difference of the original and current
* heap pointers.
*/
Meminfo
(
&
(
this
->
memSize
));
/*
* On linux systems (e.g. MCR) if built to have floating point exceptions
* turned on, invoke macro to do so
*/
#ifdef FPES_ON
unmask_std_fpes
();
#endif
this
->
paradis_ptr
=
InitHome
();
TimerInit
(
this
->
paradis_ptr
);
TimerStart
(
this
->
paradis_ptr
,
TOTAL_TIME
);
TimerStart
(
this
->
paradis_ptr
,
INITIALIZE
);
if
(
group
.
getMyRank
()
==
0
)
std
::
cout
<<
"BEGIN: Initialize ParaDiS on "
<<
nbprocs
<<
" processors."
<<
std
::
endl
;
Initialize
(
this
->
paradis_ptr
,
argc
,
const_cast
<
char
**>
(
&
argv
[
0
]),
group
.
getMPICommGroup
());
if
(
strcmp
(
this
->
paradis_ptr
->
param
->
timestepIntegrator
,
"trapezoidtimestep"
)
!=
0
)
{
LM_FATAL
(
"The timestep integrator for ParaDiS needs to be "
<<
"
\"
trapezoidtimestep
\"
in order for ParaDis to work with "
<<
"LibMultiscale. You specified
\"
timestepIntegrator = "
<<
this
->
paradis_ptr
->
param
->
timestepIntegrator
<<
"
\"
in the control file
\"
"
<<
control_filename
<<
"
\"
."
);
}
TimerStop
(
this
->
paradis_ptr
,
INITIALIZE
);
if
(
nbprocs
>
1
)
MPI_Barrier
(
this
->
paradis_ptr
->
MPI_COMM_PARADIS
);
this
->
paradis_ptr
->
param
->
timeMax
=
this
->
timestep
;
this
->
paradis_ptr
->
cycle
=
this
->
paradis_ptr
->
param
->
cycleStart
;
/*
* Perform the needed number (if any) of load-balance-only
* steps before doing the main processing loop. These steps
* perform only the minimal amount of stuff needed to
* estimate per-process load, move boundaries and migrate
* nodes among processsors to get a good initial balance.
*/
time_t
tp
;
int
initialDLBCycles
=
this
->
paradis_ptr
->
param
->
numDLBCycles
;
TimerStart
(
this
->
paradis_ptr
,
INITIALIZE
);
if
((
this
->
paradis_ptr
->
myDomain
==
0
)
&&
(
initialDLBCycles
!=
0
))
{
time
(
&
tp
);
printf
(
" +++ Beginning %d load-balancing steps at %s"
,
initialDLBCycles
,
asctime
(
localtime
(
&
tp
)));
}
while
(
this
->
paradis_ptr
->
param
->
numDLBCycles
>
0
)
{
ParadisStep
(
this
->
paradis_ptr
);
this
->
paradis_ptr
->
cycle
++
;
this
->
paradis_ptr
->
param
->
numDLBCycles
--
;
}
if
((
this
->
paradis_ptr
->
myDomain
==
0
)
&&
(
initialDLBCycles
!=
0
))
{
time
(
&
tp
);
printf
(
" +++ Completed load-balancing steps at %s"
,
asctime
(
localtime
(
&
tp
)));
}
TimerStop
(
this
->
paradis_ptr
,
INITIALIZE
);
/*
* Any time spent doing the initial DLB-only steps should
* just be attributed to initialization time, so be sure to
* reset the other timers before going into the main
* computational loop
*/
TimerInitDLBReset
(
this
->
paradis_ptr
);
/*
* The cycle number has been incremented during the initial
* load-balance steps, and previous steps so reset it to the
* proper startingvalue before entering the main processing loop.
* In paradis, cycle uses 1-based indexing.
*/
this
->
paradis_ptr
->
cycle
=
this
->
paradis_ptr
->
param
->
cycleStart
;
this
->
paradis_ptr
->
cycleLMS
=
1
;
this
->
container_dd
.
getContainerNodes
().
setGrain
(
this
->
paradis_ptr
);
// this->getContainerDisplacements().setGrain(this->paradis_ptr);
#ifndef LM_OPTIMIZED
// Debug tests.
unsigned
int
counter
=
0
;
DUMP
(
"Printing all initial ParaDiS nodes..."
,
DBG_INFO
);
for
(
auto
&&
node
:
this
->
container_dd
.
getContainerNodes
())
{
DUMP
(
"Processor "
<<
lm_my_proc_id
<<
", Node "
<<
++
counter
<<
"("
<<
node
.
position
()
<<
")"
,
DBG_INFO
);
}
#endif
// ContainerParaDiSDisplacement displs = this->getContainerDisplacements();
// DUMP("Number of ParaDiS displacement fieldpoints: "<<displs.size(),
// DBG_INFO);
// End debug tests.
if
(
this
->
paradis_ptr
->
myDomain
==
0
)
std
::
cout
<<
"END: Initialize ParaDiS on "
<<
nbprocs
<<
" processors."
<<
std
::
endl
;
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
performStep1
()
{
this
->
paradis_ptr
->
param
->
timeLMSCycle
=
0.
;
while
(
this
->
paradis_ptr
->
param
->
timeLMSCycle
<
this
->
paradis_ptr
->
param
->
timeMax
)
{
STARTTIMER
(
"ParadisStep"
);
ParadisStep
(
this
->
paradis_ptr
);
TimerClearAll
(
this
->
paradis_ptr
);
STOPTIMER
(
"ParadisStep"
);
}
this
->
paradis_ptr
->
cycleLMS
++
;
this
->
getOutput
()
->
incRelease
();
if
(
this
->
boundary
)
{
auto
c
=
GeometryManager
::
getManager
().
getGeometry
(
this
->
geom
)
->
getBoundingBox
();
auto
&&
maxval
=
c
.
getXmax
();
auto
&&
minval
=
c
.
getXmin
();
// Geometry* dd = GeometryManager::getManager().getGeometry(this->geom);
// Cube & dd_cube = dynamic_cast<Cube&>(*dd);
// to mimic strong grain boundary...
// Real xmin = this->paradis_ptr->param->xBoundMin + 10.0;
// Real xmax = this->paradis_ptr->param->xBoundMax - 10.0;
// Real ymin = this->paradis_ptr->param->yBoundMin + 10.0;
// Real ymax = this->paradis_ptr->param->yBoundMax - 10.0;
// Real zmin = this->paradis_ptr->param->zBoundMin + 10.0;
// Real zmax = this->paradis_ptr->param->zBoundMax - 10.0;
// Real xmin = this->paradis_ptr->param->xBoundMin;
// Real xmax = this->paradis_ptr->param->xBoundMax;
// Real ymin = this->paradis_ptr->param->yBoundMin;
// Real ymax = this->paradis_ptr->param->yBoundMax;
// Real zmin = this->paradis_ptr->param->zBoundMin;
// Real zmax = this->paradis_ptr->param->zBoundMax;
Real
xmin
=
minval
[
0
];
Real
xmax
=
maxval
[
0
];
Real
ymin
=
minval
[
1
];
Real
ymax
=
maxval
[
1
];
Real
zmin
=
minval
[
2
];
Real
zmax
=
maxval
[
2
];
for
(
auto
&&
node
:
this
->
container_dd
)
{
if
((
xmin
>
node
.
position
()[
0
])
||
(
node
.
position
()[
0
]
>
xmax
)
||
(
ymin
>
node
.
position
()[
1
])
||
(
node
.
position
()[
1
]
>
ymax
)
||
(
zmin
>
node
.
position
()[
2
])
||
(
node
.
position
()[
2
]
>
zmax
))
{
node
.
slaving
();
node
.
constrain
();
}
}
}
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
performStep2
()
{
// this->paradis_ptr->param->timeLMSCycle = 0.;
// while (this->paradis_ptr->param->timeLMSCycle <
// this->paradis_ptr->param->timeMax){
// ParadisStep(this->paradis_ptr);
// TimerClearAll(this->paradis_ptr);
// }
// this->paradis_ptr->cycleLMS++;
// this->dd_container.incRelease();
}
void
DomainParaDiS
::
performStep3
()
{}
void
DomainParaDiS
::
performStep4
()
{}
void
DomainParaDiS
::
performStep5
()
{}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
addSegments
(
std
::
vector
<
RefPointData
<
3
>>
&
positions
,
std
::
vector
<
RefGenericElem
<
3
>>
&
conn
,
std
::
vector
<
Real
>
&
burgers
,
std
::
vector
<
Real
>
&
normals
)
{
UInt
nb_points
=
positions
.
size
();
std
::
vector
<
Node_t
*>
created_nodes
;
std
::
vector
<
UInt
>
numNbrs
;
for
(
UInt
i
=
0
;
i
<
nb_points
;
++
i
)
{
Node_t
*
newNode
=
GetNewNativeNode
(
this
->
paradis_ptr
);
auto
&&
pos
=
positions
[
i
].
position
();
newNode
->
x
=
pos
[
0
];
newNode
->
y
=
pos
[
1
];
newNode
->
z
=
pos
[
2
];
newNode
->
native
=
1
;
newNode
->
constraint
=
7
;
newNode
->
oldvX
=
0.0
;
newNode
->
oldvY
=
0.0
;
newNode
->
oldvZ
=
0.0
;
AssignNodeToCell
(
this
->
paradis_ptr
,
newNode
);
created_nodes
.
push_back
(
newNode
);
}
numNbrs
.
resize
(
created_nodes
.
size
());
UInt
nb_elem
=
conn
.
size
();
for
(
UInt
i
=
0
;
i
<
nb_elem
;
++
i
)
{
RefGenericElem
<
3
>
&
el
=
conn
[
i
];
LM_ASSERT
(
el
.
conn
.
size
()
==
2
,
"we cannot add other things than segments to DD model"
);
for
(
UInt
n
=
0
;
n
<
2
;
++
n
)
{
UInt
nd
=
el
.
conn
[
n
];
++
numNbrs
[
nd
];
}
}
for
(
UInt
i
=
0
;
i
<
nb_points
;
++
i
)
{
Node_t
*
newNode
=
created_nodes
[
i
];
AllocNodeArms
(
newNode
,
numNbrs
[
i
]);
numNbrs
[
i
]
=
0
;
}
for
(
UInt
i
=
0
;
i
<
nb_elem
;
++
i
)
{
RefGenericElem
<
3
>
&
el
=
conn
[
i
];
UInt
nd_src
=
el
.
conn
[
0
];
UInt
nd_dst
=
el
.
conn
[
1
];
Node_t
*
srcNode
=
created_nodes
[
nd_src
];
Node_t
*
dstNode
=
created_nodes
[
nd_dst
];
UInt
&
current_arm_src
=
numNbrs
[
nd_src
];
UInt
&
current_arm_dst
=
numNbrs
[
nd_dst
];
srcNode
->
nbrTag
[
current_arm_src
].
domainID
=
dstNode
->
myTag
.
domainID
;
srcNode
->
nbrTag
[
current_arm_src
].
index
=
dstNode
->
myTag
.
index
;
srcNode
->
burgX
[
current_arm_src
]
=
burgers
[
3
*
i
+
0
];
srcNode
->
burgY
[
current_arm_src
]
=
burgers
[
3
*
i
+
1
];
srcNode
->
burgZ
[
current_arm_src
]
=
burgers
[
3
*
i
+
2
];
srcNode
->
nx
[
current_arm_src
]
=
normals
[
3
*
i
+
0
];
srcNode
->
ny
[
current_arm_src
]
=
normals
[
3
*
i
+
1
];
srcNode
->
nz
[
current_arm_src
]
=
normals
[
3
*
i
+
2
];
dstNode
->
nbrTag
[
current_arm_dst
].
domainID
=
srcNode
->
myTag
.
domainID
;
dstNode
->
nbrTag
[
current_arm_dst
].
index
=
srcNode
->
myTag
.
index
;
dstNode
->
burgX
[
current_arm_dst
]
=
-
1.
*
burgers
[
3
*
i
+
0
];
dstNode
->
burgY
[
current_arm_dst
]
=
-
1.
*
burgers
[
3
*
i
+
1
];
dstNode
->
burgZ
[
current_arm_dst
]
=
-
1.
*
burgers
[
3
*
i
+
2
];
dstNode
->
nx
[
current_arm_dst
]
=
normals
[
3
*
i
+
0
];
dstNode
->
ny
[
current_arm_dst
]
=
normals
[
3
*
i
+
1
];
dstNode
->
nz
[
current_arm_dst
]
=
normals
[
3
*
i
+
2
];
++
current_arm_src
;
++
current_arm_dst
;
}
for
(
UInt
i
=
0
;
i
<
nb_points
;
++
i
)
{
if
(
numNbrs
[
i
]
==
1
)
{
created_nodes
[
i
]
->
constraint
=
PINNED_NODE
;
// Add all pinned nodes to the list of (possible) hybrid nodes
nh_l_conInt
.
push_back
(
created_nodes
[
i
]);
}
}
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
imageForce
(
const
double
*
position
,
const
double
*
normal
,
double
*
force
)
{
/*computes the image forces on a interface given by it's normal vector and
position
needed for basic CADD
*/
unsigned
int
n_dim
=
3
;
double
totStress
[
3
][
3
];
GetFieldPointStress
(
this
->
paradis_ptr
,
position
[
0
],
position
[
1
],
position
[
2
],
totStress
);
for
(
unsigned
int
i
=
0
;
i
<
n_dim
;
++
i
)
{
force
[
i
]
=
0.
;
for
(
unsigned
int
j
=
0
;
j
<
n_dim
;
++
j
)
{
force
[
i
]
+=
totStress
[
i
][
j
]
*
normal
[
j
];
}
}
}
// CHANGE THIS?
/* -------------------------------------------------------------------------- */
/* Configure segment connecting two nodes */
void
setNodeArm
(
Node_t
*
nh_src
,
Node_t
*
nh_dst
,
const
UInt
idxArmSrc
,
const
UInt
idxArmDst
,
const
Real
*
burgersVec
,
const
Real
*
gpNormVec
)
{
/* Source node */
nh_src
->
nbrTag
[
idxArmSrc
].
domainID
=
nh_dst
->
myTag
.
domainID
;
nh_src
->
nbrTag
[
idxArmSrc
].
index
=
nh_dst
->
myTag
.
index
;
nh_src
->
burgX
[
idxArmSrc
]
=
burgersVec
[
0
];
nh_src
->
burgY
[
idxArmSrc
]
=
burgersVec
[
1
];
nh_src
->
burgZ
[
idxArmSrc
]
=
burgersVec
[
2
];
nh_src
->
nx
[
idxArmSrc
]
=
gpNormVec
[
0
];
nh_src
->
ny
[
idxArmSrc
]
=
gpNormVec
[
1
];
nh_src
->
nz
[
idxArmSrc
]
=
gpNormVec
[
2
];
/* Distributive node */
nh_dst
->
nbrTag
[
idxArmDst
].
domainID
=
nh_src
->
myTag
.
domainID
;
nh_dst
->
nbrTag
[
idxArmDst
].
index
=
nh_src
->
myTag
.
index
;
nh_dst
->
burgX
[
idxArmDst
]
=
-
burgersVec
[
0
];
nh_dst
->
burgY
[
idxArmDst
]
=
-
burgersVec
[
1
];
nh_dst
->
burgZ
[
idxArmDst
]
=
-
burgersVec
[
2
];
nh_dst
->
nx
[
idxArmDst
]
=
gpNormVec
[
0
];
nh_dst
->
ny
[
idxArmDst
]
=
gpNormVec
[
1
];
nh_dst
->
nz
[
idxArmDst
]
=
gpNormVec
[
2
];
}
/* -------------------------------------------------------------------------- */
/* Orthogonal projection of a node onto the glide plane */
inline
void
projectNodeOnGP
(
Node_t
*
nh_src
,
// source node
Node_t
*
nh_dst
,
// inheriting node
const
UInt
idxArmSrc
)
// arm index of source node
// connecting to the
// inheriting node
{
Real
alpha
=
nh_src
->
nx
[
idxArmSrc
]
*
(
nh_src
->
x
-
nh_dst
->
x
)
+
nh_src
->
ny
[
idxArmSrc
]
*
(
nh_src
->
y
-
nh_dst
->
y
)
+
nh_src
->
nz
[
idxArmSrc
]
*
(
nh_src
->
z
-
nh_dst
->
z
);
nh_dst
->
x
+=
alpha
*
nh_src
->
nx
[
idxArmSrc
];
nh_dst
->
y
+=
alpha
*
nh_src
->
ny
[
idxArmSrc
];
nh_dst
->
z
+=
alpha
*
nh_src
->
nz
[
idxArmSrc
];
}
/* -------------------------------------------------------------------------- */
/* Add hybrid segments to an existing DD line crossing an interface */
// NOTE: test version --> the method assumes that the mesh struct carries only 1
// dislocation!!
// NOTE: it is further assumed that the node of hybrid DD line in the continuum
// connecting to a node in the atomistic domain remains the nearest node
// to the interface during a simulation!
// NOTE: there is no particular check performed if "continuum" DD nodes
// previously connected to a DD node in the atomistic domain retain the
// this connectivity!!!
// NOTE:
void
DomainParaDiS
::
addHybSegments
(
std
::
vector
<
RefPointData
<
3
>>
&
positions
,
// DD nodes in an atomistic domain
std
::
vector
<
RefGenericElem
<
3
>>
&
conn
,
// connectivity between DD nodes in the atomistic domain
Real
threshold
)
// max. distance between two adjacent DD nodes possibly
// generating a hybrid segment
{
UInt
nb_points
=
positions
.
size
();
Real
dist
;
Real
pos_ato
[
3
],
pos_con
[
3
];
Real
burgersVec
[
3
],
gpNormVec
[
3
];
// Containers/Lists
std
::
vector
<
UInt
>
numNbrs
;
// ParaDis
Node_t
*
nh
,
*
nh_conInt1
=
NULL
,
*
nh_conInt2
=
NULL
;
// ParaDis - Iterators
std
::
vector
<
Node_t
*>::
iterator
nh_it
;
/* ------------------------------------------------------------------------ */
// Delete all existing DD segments in the atomistic domain
// NOTE: the FreeNode() function doesn't remove the segments associated with
// the node to delete
// --> more precisely, they are not deallocated
// --> this is of no particular issue since one can call
// ReallocNodeArms() in step 3)
for
(
nh_it
=
created_nodes
.
begin
();
nh_it
!=
created_nodes
.
end
();
++
nh_it
)
{
FreeNode
(
paradis_ptr
,
(
*
nh_it
)
->
myTag
.
index
);
}
created_nodes
.
clear
();
/* ------------------------------------------------------------------------ */
/* 1) Import new nodes in ParaDis */
for
(
UInt
i
=
0
;
i
<
nb_points
;
++
i
)
{
nh
=
GetNewNativeNode
(
this
->
paradis_ptr
);
auto
&&
pos
=
positions
[
i
].
position
();
nh
->
x
=
pos
[
0
];
nh
->
y
=
pos
[
1
];
nh
->
z
=
pos
[
2
];
nh
->
native
=
1
;
nh
->
constraint
=
PINNED_NODE
;
// nodes in an atomistic domain remain always fixed
nh
->
oldvX
=
0.0
;
nh
->
oldvY
=
0.0
;
nh
->
oldvZ
=
0.0
;
AssignNodeToCell
(
this
->
paradis_ptr
,
nh
);
created_nodes
.
push_back
(
nh
);
}
numNbrs
.
resize
(
created_nodes
.
size
());
/* ------------------------------------------------------------------------ */
/* 2) Find nearest neighbor(s) in the continuum */
// 2.1) First node
nh
=
created_nodes
[
0
];
pos_ato
[
0
]
=
nh
->
x
;
pos_ato
[
1
]
=
nh
->
y
;
pos_ato
[
2
]
=
nh
->
z
;
for
(
nh_it
=
this
->
nh_l_conInt
.
begin
();
nh_it
!=
this
->
nh_l_conInt
.
end
();
++
nh_it
)
{
// Position of "continuum" node
pos_con
[
0
]
=
(
*
nh_it
)
->
x
;
pos_con
[
1
]
=
(
*
nh_it
)
->
y
;
pos_con
[
2
]
=
(
*
nh_it
)
->
z
;
// Compute the distance between the continuum and the
dist
=
sqrt
((
pos_ato
[
0
]
-
pos_con
[
0
])
*
(
pos_ato
[
0
]
-
pos_con
[
0
])
+
(
pos_ato
[
1
]
-
pos_con
[
1
])
*
(
pos_ato
[
1
]
-
pos_con
[
1
])
+
(
pos_ato
[
2
]
-
pos_con
[
2
])
*
(
pos_ato
[
2
]
-
pos_con
[
2
]));
// If distance is smaller than a prescribed threshold
// --> mark the node as a neighbor of the first "atomistic" DD node
if
(
dist
<
threshold
)
{
nh_conInt1
=
*
nh_it
;
numNbrs
[
0
]
=
1
;
break
;
}
}
// 2.2) Last node
nh
=
created_nodes
[
nb_points
-
1
];
pos_ato
[
0
]
=
nh
->
x
;
pos_ato
[
1
]
=
nh
->
y
;
pos_ato
[
2
]
=
nh
->
z
;
for
(
nh_it
=
this
->
nh_l_conInt
.
begin
();
nh_it
!=
this
->
nh_l_conInt
.
end
();
++
nh_it
)
{
// Position of "continuum" node
pos_con
[
0
]
=
(
*
nh_it
)
->
x
;
pos_con
[
1
]
=
(
*
nh_it
)
->
y
;
pos_con
[
2
]
=
(
*
nh_it
)
->
z
;
// Compute the distance between the continuum and the
dist
=
sqrt
((
pos_ato
[
0
]
-
pos_con
[
0
])
*
(
pos_ato
[
0
]
-
pos_con
[
0
])
+
(
pos_ato
[
1
]
-
pos_con
[
1
])
*
(
pos_ato
[
1
]
-
pos_con
[
1
])
+
(
pos_ato
[
2
]
-
pos_con
[
2
])
*
(
pos_ato
[
2
]
-
pos_con
[
2
]));
// If distance is smaller than a prescribed threshold
// --> mark the node as a neighbor of the last "atomistic" DD node
if
(
dist
<
threshold
)
{
nh_conInt2
=
*
nh_it
;
numNbrs
[
nb_points
-
1
]
=
1
;
break
;
}
}
/* ------------------------------------------------------------------------ */
/* 3) Import connectivity and allocate memory for segment info */
UInt
nb_elem
=
conn
.
size
();
for
(
UInt
i
=
0
;
i
<
nb_elem
;
++
i
)
{
RefGenericElem
<
3
>
&
el
=
conn
[
i
];
LM_ASSERT
(
el
.
conn
.
size
()
==
2
,
"we cannot add other things than segments to DD model"
);
for
(
UInt
n
=
0
;
n
<
2
;
++
n
)
{
UInt
nd
=
el
.
conn
[
n
];
++
numNbrs
[
nd
];
}
}
for
(
UInt
i
=
0
;
i
<
nb_points
;
++
i
)
{
ReallocNodeArms
(
created_nodes
[
i
],
numNbrs
[
i
]);
// AllocNodeArms(newNode, numNbrs[i]);
numNbrs
[
i
]
=
0
;
}
/* ------------------------------------------------------------------------ */
/* 4) Connect atomistic and continuum DD line */
// 4.1) First interface node in the continuum connects to the atomistic DD
// line
if
((
nh_conInt1
!=
NULL
)
&&
(
nh_conInt2
==
NULL
))
{
// Get the Burgers vector of the dislocation
burgersVec
[
0
]
=
nh_conInt1
->
burgX
[
0
];
burgersVec
[
1
]
=
nh_conInt1
->
burgY
[
0
];
burgersVec
[
2
]
=
nh_conInt1
->
burgZ
[
0
];
// Get the normal vector of the glide plane
gpNormVec
[
0
]
=
nh_conInt1
->
nx
[
0
];
gpNormVec
[
1
]
=
nh_conInt1
->
ny
[
0
];
gpNormVec
[
2
]
=
nh_conInt1
->
nz
[
0
];
// Re-allocate memory for the hybrid segment if necessary (--> should only
// be done once before the first iteration)
if
(
nh_conInt1
->
numNbrs
<
2
)
{
ReallocNodeArms
(
nh_conInt1
,
2
);
nh_conInt1
->
constraint
=
UNCONSTRAINED
;
}
// Establish connectivity between nodes in the continuum and the atomistic
// domain ...
setNodeArm
(
nh_conInt1
,
created_nodes
[
0
],
1
,
0
,
burgersVec
,
gpNormVec
);
// // ... and set homogeneous Dirichlet conditions on the other node
// (--> u=0 <=> pinned node)
// created_nodes[nb_points-1]->constraint = PINNED_NODE;
// Project atomistic DD node on the glide plane
projectNodeOnGP
(
nh_conInt1
,
created_nodes
[
0
],
1
);
++
numNbrs
[
0
];
}
// 4.2) Last interface node in the continuum connects to the atomistic DD line
else
if
((
nh_conInt1
==
NULL
)
&&
(
nh_conInt2
!=
NULL
))
{
// Get the Burgers vector of the dislocation
burgersVec
[
0
]
=
nh_conInt2
->
burgX
[
0
];
burgersVec
[
1
]
=
nh_conInt2
->
burgY
[
0
];
burgersVec
[
2
]
=
nh_conInt2
->
burgZ
[
0
];
// Get the normal vector of the glide plane
gpNormVec
[
0
]
=
nh_conInt2
->
nx
[
0
];
gpNormVec
[
1
]
=
nh_conInt2
->
ny
[
0
];
gpNormVec
[
2
]
=
nh_conInt2
->
nz
[
0
];
// Re-allocate memory for the hybrid segment if necessary (--> should only
// be done once before the first iteration)
if
(
nh_conInt2
->
numNbrs
<
2
)
{
ReallocNodeArms
(
nh_conInt2
,
2
);
nh_conInt2
->
constraint
=
UNCONSTRAINED
;
}
// Establish connectivity between nodes in the continuum and the atomistic
// domain ...
setNodeArm
(
nh_conInt2
,
created_nodes
[
nb_points
-
1
],
1
,
0
,
burgersVec
,
gpNormVec
);
// // ... and set homogeneous Dirichlet conditions on the other node
// (--> u=0 <=> pinned node)
// created_nodes[0]->constraint = PINNED_NODE;
// Project atomistic DD node on the glide plane
projectNodeOnGP
(
nh_conInt2
,
created_nodes
[
nb_points
-
1
],
1
);
++
numNbrs
[
nb_points
-
1
];
}
// 4.3) Both interfaces nodes in the continuum connects to the atomistic DD
// line
else
if
((
nh_conInt1
!=
NULL
)
&&
(
nh_conInt2
!=
NULL
))
{
// Get the Burgers vector of the dislocation
burgersVec
[
0
]
=
nh_conInt1
->
burgX
[
0
];
burgersVec
[
1
]
=
nh_conInt1
->
burgY
[
0
];
burgersVec
[
2
]
=
nh_conInt1
->
burgZ
[
0
];
// Get the normal vector of the glide plane
gpNormVec
[
0
]
=
nh_conInt1
->
nx
[
0
];
gpNormVec
[
1
]
=
nh_conInt1
->
ny
[
0
];
gpNormVec
[
2
]
=
nh_conInt1
->
nz
[
0
];
// Re-allocate memory for the hybrid segment if necessary (--> should only
// be done once before the first iteration)
if
(
nh_conInt1
->
numNbrs
<
2
)
{
ReallocNodeArms
(
nh_conInt1
,
2
);
nh_conInt1
->
constraint
=
UNCONSTRAINED
;
}
if
(
nh_conInt2
->
numNbrs
<
2
)
{
ReallocNodeArms
(
nh_conInt2
,
2
);
nh_conInt2
->
constraint
=
UNCONSTRAINED
;
}
// Establish connectivity between nodes in the continuum and the atomistic
// domain
setNodeArm
(
nh_conInt1
,
created_nodes
[
0
],
1
,
0
,
burgersVec
,
gpNormVec
);
setNodeArm
(
nh_conInt2
,
created_nodes
[
nb_points
-
1
],
1
,
0
,
burgersVec
,
gpNormVec
);
// Project atomistic DD nodes on the glide plane
projectNodeOnGP
(
nh_conInt1
,
created_nodes
[
0
],
1
);
projectNodeOnGP
(
nh_conInt1
,
created_nodes
[
nb_points
-
1
],
1
);
++
numNbrs
[
0
];
++
numNbrs
[
nb_points
-
1
];
}
else
{
LM_FATAL
(
std
::
endl
<<
"Hybrid dislocation failed to connect!!"
<<
std
::
endl
);
}
/* ------------------------------------------------------------------------ */
/* 5) Now connect the remainder of the hybrid dislocation line in the
* atomistic domain */
for
(
UInt
i
=
0
;
i
<
nb_elem
;
++
i
)
{
// Get connectivity
RefGenericElem
<
3
>
&
el
=
conn
[
i
];
UInt
idxNhSrc
=
el
.
conn
[
0
];
UInt
idxNhDst
=
el
.
conn
[
1
];
UInt
&
idxArmSrc
=
numNbrs
[
idxNhSrc
];
UInt
&
idxArmDst
=
numNbrs
[
idxNhDst
];
// Update arm information for each node of segment i
setNodeArm
(
created_nodes
[
idxNhSrc
],
created_nodes
[
idxNhDst
],
idxArmSrc
,
idxArmDst
,
burgersVec
,
gpNormVec
);
// Project atomistic DD node on the glide plane
projectNodeOnGP
(
created_nodes
[
idxNhSrc
],
created_nodes
[
idxNhDst
],
idxArmSrc
);
++
idxArmSrc
;
++
idxArmDst
;
}
}
/* -------------------------------------------------------------------------- */
UInt
DomainParaDiS
::
getNBNodes
()
{
return
this
->
container_dd
.
getContainerNodes
().
size
();
}
/* -------------------------------------------------------------------------- */
Real
DomainParaDiS
::
getShearModulus
()
{
return
this
->
paradis_ptr
->
param
->
shearModulus
;
}
/* -------------------------------------------------------------------------- */
VectorProxy
<
3u
>
DomainParaDiS
::
getXmin
()
{
return
VectorProxy
<
3u
>
(
this
->
paradis_ptr
->
param
->
xBoundMin
,
this
->
paradis_ptr
->
param
->
yBoundMin
,
this
->
paradis_ptr
->
param
->
zBoundMin
);
}
/* -------------------------------------------------------------------------- */
VectorProxy
<
3u
>
DomainParaDiS
::
getXmax
()
{
return
VectorProxy
<
3u
>
(
this
->
paradis_ptr
->
param
->
xBoundMax
,
this
->
paradis_ptr
->
param
->
yBoundMax
,
this
->
paradis_ptr
->
param
->
zBoundMax
);
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
freenodearm
(
UInt
index
)
{
Node_t
node
=
*
paradis_ptr
->
nodeKeys
[
index
];
FreeNodeArms
(
&
node
);
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
freenode
(
UInt
index
)
{
Node_t
node
=
*
paradis_ptr
->
nodeKeys
[
index
];
FreeNodeArms
(
&
node
);
}
/* -------------------------------------------------------------------------- */
bool
DomainParaDiS
::
is_it_bound
()
{
// ContainerParaDiSNode::iterator dbgNodes =
// this->getContainerNodes().begin();
// UInt counter = 0;
// for(RefParaDiSNode node=dbgNodes.getFirst(); !dbgNodes.end();
// node=dbgNodes.getNext()){
// ++counter;
// }
return
this
->
boundary
;
}
// CHANGE THIS?
/* -------------------------------------------------------------------------- */
Cube
DomainParaDiS
::
getDomainBoundingBox
()
{
Cube
c
;
c
.
getXmin
()
=
VectorView
<
3
>
(
this
->
paradis_ptr
->
param
->
minCoordinates
);
c
.
getXmax
()
=
VectorView
<
3
>
(
this
->
paradis_ptr
->
param
->
maxCoordinates
);
return
c
;
}
/* -------------------------------------------------------------------------- */
Cube
DomainParaDiS
::
getSubDomainBoundingBox
()
{
LM_TOIMPLEMENT
;
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
setDomainBoundingBox
(
const
Cube
&
)
{
LM_TOIMPLEMENT
;
}
/* -------------------------------------------------------------------------- */
void
DomainParaDiS
::
setSubDomainBoundingBox
(
const
Cube
&
)
{
LM_TOIMPLEMENT
;
}
// void getDomainDimensions(Real *xmin, Real *xmax) {
// };
/* -------------------------------------------------------------------------- */
/* LMDESC PARADIS
This section is used to configure the ParaDiS plugin
*/
/* LMEXAMPLE
Section MultiScale AtomsUnits\\
...\\
GEOMETRY 1 CUBE BBOX -10 10 -10 10 -10 10
MODEL PARADIS dislomodel \\
...\\
endSection\\ \\
Section PARADIS:dislomodel RealUnits \\
NUMDLBSTEPS 0 \\
CONTROLFILE frank_read_src.ctrl \\
DATAFILE frank_read_src.data \\
TIMESTEP 1e-10 \\
DOMAIN_GEOMETRY 1 \\
endSection
*/
/* LMHERITATE domain_dd */
void
DomainParaDiS
::
declareParams
()
{
DomainDD
<
ContainerParaDiSNode
,
3
>::
declareParams
();
/* LMKEYWORD TIMESTEP
Set the maximum time step used for ParaDiS
*/
this
->
parseKeyword
(
"TIMESTEP"
,
timestep
);
/* LMKEYWORD NUMDLBSTEPS
Number of dynamic load balancing steps
*/
this
->
parseKeyword
(
"NUMDLBSTEPS"
,
numDLBCycles
);
/* LMKEYWORD CONTROLFILE
Path to the ParaDiS control file
*/
this
->
parseKeyword
(
"CONTROLFILE"
,
control_filename
);
/* LMKEYWORD DATAFILE
Path to the ParaDis data file
*/
this
->
parseKeyword
(
"DATAFILE"
,
data_filename
);
/* LMKEYWORD MOBILITYFILE
Mobility input acquired from MD
*/
this
->
parseKeyword
(
"MOBILITYFILE"
,
mobility_filename
,
""
);
/* LMKEYWORD BOUNDARY
Constraint boundary to mimic grain boundary
*/
this
->
parseTag
(
"BOUNDARY"
,
this
->
boundary
,
false
);
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
Event Timeline
Log In to Comment