Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81138049
ATC_Transfer.cpp
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 4, 09:32
Size
46 KB
Mime Type
text/x-c
Expires
Fri, Sep 6, 09:32 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
20393886
Attached To
rLAMMPS lammps
ATC_Transfer.cpp
View Options
// ATC_Transfer headers
#include "ATC_Transfer.h"
#include "FE_Engine.h"
#include "Array.h"
#include "Array2D.h"
#include "ATC_Error.h"
#include "CG.h"
#include "XT_Function.h"
#include "PrescribedDataManager.h"
#include "TimeIntegrator.h"
#include "PhysicsModel.h"
#include "PhysicsModelThermal.h"
#include "PhysicsModelTwoTemperature.h"
namespace
ATC
{
ATC_Transfer
::
ATC_Transfer
(
void
)
:
lammpsInterface_
(
LammpsInterface
::
instance
()),
physicsModel_
(
NULL
),
feEngine_
(
NULL
),
prescribedDataMgr_
(
NULL
),
simTime_
(
0.0
),
stepCounter_
(
0
),
sampleCounter_
(
0
),
outputFrequency_
(
0
),
sampleFrequency_
(
1
),
outputFrequencyAtom_
(
0
),
nLocal_
(
0
),
nLocalTotal_
(
0
),
nLocalGhost_
(
0
),
nLocalMask_
(
0
),
nInternal_
(
0
),
nGhost_
(
0
),
nLocalLambda_
(
0
),
equilibriumStart_
(
false
),
mdOutputManager_
(),
bndyIntType_
(
NO_QUADRATURE
),
bndyFaceSet_
(
NULL
),
globalSwitch_
(
0
),
resetSwitch_
(
false
),
atomToElementMapType_
(
LAGRANGIAN
),
atomToElementMapFrequency_
(
0
),
neighborResetFrequency_
(
0
),
scalarFlag_
(
0
),
vectorFlag_
(
0
),
sizeVector_
(
0
),
globalFreq_
(
0
),
extScalar_
(
0
),
extVector_
(
0
),
extList_
(
NULL
),
timeFilterManager_
(
this
),
groupbit_
(
0
),
groupbitGhost_
(
0
),
atomSwitch_
(
false
),
initialized_
(
false
),
atomQuadForInternal_
(
true
),
useLocalizedLambda_
(
false
),
useLumpedLambda_
(
false
),
regionID_
(
-
1
),
atomWeightType_
(
LATTICE
),
xref_
(
NULL
),
readXref_
(
false
),
timeIntegrator_
(
NULL
),
extrinsicModelManager_
(
this
),
useRestart_
(
false
),
trackCharge_
(
false
)
{
// Defaults
grow_arrays
(
lammpsInterface_
->
nmax
());
fieldMask_
.
reset
(
NUM_FIELDS
,
NUM_FLUX
);
fieldMask_
=
false
;
// fe_engine
feEngine_
=
new
FE_Engine
(
this
);
// check to see if lammps has any charges
double
*
atomicCharge
=
lammpsInterface_
->
atom_charge
();
if
(
atomicCharge
)
trackCharge_
=
true
;
}
ATC_Transfer
::~
ATC_Transfer
()
{
lammpsInterface_
->
destroy_2d_double_array
(
xref_
);
if
(
feEngine_
)
delete
feEngine_
;
if
(
physicsModel_
)
delete
physicsModel_
;
if
(
prescribedDataMgr_
)
delete
prescribedDataMgr_
;
}
//--------------------------------------------------
// pack_fields
// bundle all allocated field matrices into a list
// for output needs
//--------------------------------------------------
void
ATC_Transfer
::
pack_fields
(
OUTPUT_LIST
&
data
)
{
map
<
FieldName
,
int
>::
const_iterator
field
;
for
(
field
=
fieldSizes_
.
begin
();
field
!=
fieldSizes_
.
end
();
field
++
)
{
FieldName
thisField
=
field
->
first
;
string
fieldName
=
field_to_string
(
thisField
);
string
matrixName
;
// copy all fields from ATC_Transfer.h
matrixName
=
"fields_"
+
fieldName
;
data
[
matrixName
]
=
&
fields_
[
thisField
];
matrixName
=
"dot_fields_"
+
fieldName
;
data
[
matrixName
]
=
&
dot_fields_
[
thisField
];
matrixName
=
"ddot_fields_"
+
fieldName
;
data
[
matrixName
]
=
&
ddot_fields_
[
thisField
];
matrixName
=
"dddot_fields_"
+
fieldName
;
data
[
matrixName
]
=
&
dddot_fields_
[
thisField
];
matrixName
=
"dot_fieldsMD_"
+
fieldName
;
data
[
matrixName
]
=
&
dot_fieldsMD_
[
thisField
];
matrixName
=
"ddot_fieldsMD_"
+
fieldName
;
data
[
matrixName
]
=
&
ddot_fieldsMD_
[
thisField
];
matrixName
=
"dot_dot_fieldsMD_"
;
data
[
matrixName
]
=
&
dot_dot_fieldsMD_
[
thisField
];
matrixName
=
"fieldNdOld_"
+
fieldName
;
data
[
matrixName
]
=
&
fieldNdOld_
[
thisField
];
matrixName
=
"fieldNdFiltered_"
+
fieldName
;
data
[
matrixName
]
=
&
fieldNdFiltered_
[
thisField
];
matrixName
=
"fieldRateNdOld_"
;
data
[
matrixName
]
=
&
fieldRateNdOld_
[
thisField
];
matrixName
=
"fieldRateNdFiltered_"
+
fieldName
;
data
[
matrixName
]
=
&
fieldRateNdFiltered_
[
thisField
];
matrixName
=
"dot_fieldRateNdOld_"
+
fieldName
;
data
[
matrixName
]
=
&
dot_fieldRateNdOld_
[
thisField
];
matrixName
=
"dot_fieldRateNdFiltered_"
+
fieldName
;
data
[
matrixName
]
=
&
dot_fieldRateNdFiltered_
[
thisField
];
matrixName
=
"auxStorage_"
+
fieldName
;
data
[
matrixName
]
=
&
auxStorage_
[
thisField
];
}
}
//--------------------------------------------------
// write_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void
ATC_Transfer
::
write_restart_data
(
string
fileName
,
OUTPUT_LIST
&
data
)
{
pack_fields
(
data
);
feEngine_
->
write_restart_file
(
fileName
,
&
data
);
}
//--------------------------------------------------
// read_restart_file
// bundle matrices that need to be saved and call
// fe_engine to write the file
//--------------------------------------------------
void
ATC_Transfer
::
read_restart_data
(
string
fileName
,
OUTPUT_LIST
&
data
)
{
pack_fields
(
data
);
feEngine_
->
read_restart_file
(
fileName
,
&
data
);
}
//--------------------------------------------------
// Interactions with LAMMPS fix commands
// parse input command and pass on to finite element engine
// or physics specific transfers if necessary
// revert to physics-specific transfer if no command matches input
// first keyword is unique to particular class
// base class keyword matching must apply to ALL physics
// order: derived, base, owned objects
//--------------------------------------------------
bool
ATC_Transfer
::
modify
(
int
narg
,
char
**
arg
)
{
FieldName
thisField
;
int
thisIndex
;
int
argIdx
;
bool
match
=
false
;
if
(
strcmp
(
arg
[
0
],
"transfer"
)
==
0
)
{
/*! \page man_transfer_output fix_modify AtC transfer output
\section syntax
fix_modify AtC transfer output <filename_prefix> <frequency>
[text | vector_components | tensor_components ]
- filename_prefix (string) = prefix for data files
- frequency (integer) = frequency of output in time-steps
- options (keyword/s): \n
text = creates text output as well as binary Ensight output \n
vector_components = outputs vectors as scalar components \n
tensor_components = outputs tensor as scalar components
(use this for Paraview)\n
\section examples
<TT> fix_modify AtC transfer output heatFE 100 </TT> \n
<TT> fix_modify AtC transfer output hardyFE 1 text tensor_components </TT> \n
\section description
Creates (binary, "gold" format) Ensight output of nodal/mesh data
which is transfer/physics specific.
\section restrictions
\section related
see \ref man_fix_atc
\section default
none
*/
if
(
strcmp
(
arg
[
1
],
"output"
)
==
0
)
{
if
(
strcmp
(
arg
[
2
],
"nodeset"
)
==
0
)
{
string
nset
=
arg
[
3
];
if
(
strcmp
(
arg
[
4
],
"sum"
)
==
0
)
{
string
field
=
arg
[
5
];
pair
<
string
,
FieldName
>
id
(
nset
,
string_to_field
(
field
));
nsetCompute_
[
id
]
=
0.0
;
match
=
true
;
}
}
else
{
outputPrefix_
=
arg
[
2
];
outputFrequency_
=
atoi
(
arg
[
3
]);
bool
text_output
=
false
,
vect_comp
=
false
,
tensor_comp
=
false
;
int
rank
=
lammpsInterface_
->
comm_rank
();
for
(
int
i
=
4
;
i
<
narg
;
++
i
)
{
if
(
strcmp
(
arg
[
i
],
"text"
)
==
0
)
text_output
=
true
;
else
if
(
strcmp
(
arg
[
i
],
"vector_components"
)
==
0
)
vect_comp
=
true
;
else
if
(
strcmp
(
arg
[
i
],
"tensor_components"
)
==
0
)
tensor_comp
=
true
;
}
if
(
outputFrequency_
>
0
)
{
if
(
text_output
)
{
feEngine_
->
initialize_output
(
rank
,
outputPrefix_
,
GNUPLOT
);
if
(
rank
==
0
)
cout
<<
" ATC:: Warning : text output can create _LARGE_ files
\n
"
;
}
else
feEngine_
->
initialize_output
(
rank
,
outputPrefix_
);
if
(
vect_comp
)
feEngine_
->
output_manager
()
->
set_option
(
OUTPUT_VECTOR_COMPONENTS
,
true
);
if
(
tensor_comp
)
feEngine_
->
output_manager
()
->
set_option
(
OUTPUT_TENSOR_COMPONENTS
,
true
);
}
match
=
true
;
}
}
/*! \page man_transfer_atomic_output fix_modify AtC transfer atomic_output
\section syntax
fix_modify AtC transfer atomic_output <filename_prefix> <frequency>
[text]
- filename_prefix (string) = prefix for data files
- frequency (integer) = frequency of output in time-steps
- option "text" = creates a text version of the output as well as
the binary Ensight output
\section examples
<TT> fix_modify AtC transfer atomic_output heatMD 100 </TT> \n
<TT> fix_modify AtC transfer atomic_output nanoMD 1 text </TT> \n
\section description
Creates (binary, "gold" format) Ensight output of atom (point) data
which is transfer/physics specific.
\section restrictions
\section related
\section default
none
*/
else
if
(
strcmp
(
arg
[
1
],
"atomic_output"
)
==
0
)
{
outputPrefixAtom_
=
arg
[
2
];
outputFrequencyAtom_
=
atoi
(
arg
[
3
]);
if
(
outputFrequencyAtom_
>
0
)
{
if
(
narg
==
5
&&
strcmp
(
arg
[
4
],
"text"
)
==
0
)
{
mdOutputManager_
.
initialize
(
outputPrefixAtom_
,
GNUPLOT
);
cout
<<
" ATC:: Warning : text output can create _LARGE_ files
\n
"
;
}
else
mdOutputManager_
.
initialize
(
outputPrefixAtom_
);
}
match
=
true
;
}
/*! \page man_transfer_internal fix_modify AtC transfer internal
\section syntax
fix_modify AtC transfer internal type <atom-type-id>
- <atom-type-id> = type id for non-fixed atoms internal to the FE mesh
\section examples
<TT> fix_modify AtC transfer internal type internal_atoms </TT>
\section description
Command to define the atoms to couple to finite elements.
This definition is required for a AtC simulation.
\section restrictions
\section related
see \ref man_transfer_boundary
\section default
none
*/
else
if
((
strcmp
(
arg
[
1
],
"internal"
)
==
0
)
&&
(
strcmp
(
arg
[
2
],
"type"
)
==
0
))
{
int
igroup
=
lammpsInterface_
->
find_group
(
arg
[
3
]);
groupbit_
|=
lammpsInterface_
->
group_bit
(
igroup
);
igroups_
.
insert
(
igroup
);
match
=
true
;
}
/*! \page man_transfer_boundary fix_modify AtC transfer boundary
\section syntax
fix_modify AtC transfer boundary type <atom-type-id>
- <atom-type-id> = type id for atoms that represent a ficticious
boundary internal to the FE mesh
\section examples
<TT> fix_modify AtC transfer boundary type ghost_atoms </TT>
\section description
Command to define the atoms that represent the ficticious
boundary internal to the FE mesh. For fully overlapped MD/FE
domains with periodic boundary conditions no boundary atoms should
be defined.
\section restrictions
\section related
see \ref man_transfer_internal
\section default
none
*/
else
if
((
strcmp
(
arg
[
1
],
"boundary"
)
==
0
)
&&
(
strcmp
(
arg
[
2
],
"type"
)
==
0
))
{
int
igroup
=
lammpsInterface_
->
find_group
(
arg
[
3
]);
groupbitGhost_
|=
lammpsInterface_
->
group_bit
(
igroup
);
igroupsGhost_
.
insert
(
igroup
);
match
=
true
;
}
/*! \page man_internal_quadrature fix_modify AtC transfer internal_quadrature
\section syntax
fix_modify atc transfer internal_quadrature < on | off >
\section examples
<TT> fix_modify atc transfer internal_quadrature off </TT>
\section description
Command use or not use atomic quadrature on internal elements
fully filled with atoms. By turning the internal quadrature off
these elements do not contribute to the governing PDE and the fields
at the internal nodes follow the weighted averages of the atomic data.
\section restrictions
\section related
\section default
on
*/
else
if
(
strcmp
(
arg
[
1
],
"internal_quadrature"
)
==
0
)
{
if
(
strcmp
(
arg
[
2
],
"on"
)
==
0
)
{
atomQuadForInternal_
=
true
;
match
=
true
;
}
else
if
(
strcmp
(
arg
[
2
],
"off"
)
==
0
)
{
atomQuadForInternal_
=
false
;
regionID_
=
-
1
;
match
=
true
;
}
else
if
(
strcmp
(
arg
[
3
],
"off"
)
==
0
)
{
for
(
regionID_
=
0
;
regionID_
<
lammpsInterface_
->
nregion
();
regionID_
++
)
if
(
strcmp
(
arg
[
2
],
lammpsInterface_
->
region_name
(
regionID_
))
==
0
)
break
;
if
(
regionID_
<
lammpsInterface_
->
nregion
())
{
atomQuadForInternal_
=
false
;
match
=
true
;
}
else
regionID_
=
-
1
;
}
}
/*! \page man_initial fix_modify AtC transfer initial
\section syntax
fix_modify AtC transfer initial <field> <nodeset> <constant | function>
- <field> = field name valid for type of physics, temperature | electron_temperature
- <nodeset> = name of set of nodes to apply boundary condition
- <constant | function> = value or name of function followed by its
parameters
\section examples
<TT> fix_modify atc transfer initial temperature groupNAME 10. </TT>
\section description
Sets the initial values for the specified field at the specified nodes.
\section restrictions
keyword 'all' reserved in nodeset name
\section related
see \ref man_transfer_internal
\section default
none
*/
// set initial conditions
else
if
(
strcmp
(
arg
[
1
],
"initial"
)
==
0
)
{
argIdx
=
2
;
get_field
(
arg
,
argIdx
,
thisField
,
thisIndex
);
string
nsetName
(
arg
[
argIdx
++
]);
XT_Function
*
f
=
NULL
;
// parse constant
if
(
narg
==
argIdx
+
1
)
{
f
=
XT_Function_Mgr
::
instance
()
->
get_constant_function
(
atof
(
arg
[
argIdx
]));
}
// parse function
else
{
f
=
XT_Function_Mgr
::
instance
()
->
get_function
(
&
(
arg
[
argIdx
]),
narg
-
argIdx
);
}
prescribedDataMgr_
->
fix_initial_field
(
nsetName
,
thisField
,
thisIndex
,
f
);
match
=
true
;
}
/*! \page man_fix_nodes fix_modify AtC transfer fix
\section syntax
fix_modify AtC transfer fix <field> <nodeset> <constant | function>
- <field> = field name valid for type of physics
- <nodeset> = name of set of nodes to apply boundary condition
- <constant | function> = value or name of function followed by its
parameters
\section examples
<TT> fix_modify AtC transfer fix temperature groupNAME 10. </TT> \n
<TT> fix_modify AtC transfer fix temperature groupNAME 0 0 0 10.0 0 0 1.0 </TT> \n
\section description
Creates a constraint on the values of the specified field at specified nodes.
\section restrictions
keyword 'all' reserved in nodeset name
\section related
see \ref man_unfix_nodes
\section default
none
*/
// fix and unfix nodes
else
if
(
strcmp
(
arg
[
1
],
"fix"
)
==
0
)
{
argIdx
=
2
;
get_field
(
arg
,
argIdx
,
thisField
,
thisIndex
);
string
nsetName
(
arg
[
argIdx
++
]);
XT_Function
*
f
=
NULL
;
// parse constant
if
(
narg
==
argIdx
+
1
)
{
f
=
XT_Function_Mgr
::
instance
()
->
get_constant_function
(
atof
(
arg
[
argIdx
]));
}
// parse function
else
{
f
=
XT_Function_Mgr
::
instance
()
->
get_function
(
&
(
arg
[
argIdx
]),
narg
-
argIdx
);
}
prescribedDataMgr_
->
fix_field
(
nsetName
,
thisField
,
thisIndex
,
f
);
match
=
true
;
}
/*! \page man_unfix_nodes fix_modify AtC transfer unfix
\section syntax
fix_modify AtC transfer unfix <field> <nodeset>
- <field> = field name valid for type of physics
- <nodeset> = name of set of nodes
\section examples
<TT> fix_modify AtC transfer unfix temperature groupNAME </TT>
\section description
Removes constraint on field values for specified nodes.
\section restrictions
keyword 'all' reserved in nodeset name
\section related
see \ref man_fix_nodes
\section default
none
*/
else
if
(
strcmp
(
arg
[
1
],
"unfix"
)
==
0
)
{
argIdx
=
2
;
get_field
(
arg
,
argIdx
,
thisField
,
thisIndex
);
string
nsetName
(
arg
[
argIdx
++
]);
prescribedDataMgr_
->
unfix_field
(
nsetName
,
thisField
,
thisIndex
);
match
=
true
;
}
/*! \page man_source fix_modify AtC transfer source
\section syntax
fix_modify AtC transfer source <field> <element_set> <value | function>
- <field> = field name valid for type of physics
- <element_set> = name of set of elements
\section examples
<TT> fix_modify atc transfer source temperature middle temporal_ramp 10. 0. </TT>
\section description
Add domain sources to the mesh. The units are consistent with LAMMPS's
units for mass, length and time and are defined by the PDE being solved,
e.g. for thermal transfer the balance equation is for energy and source
is energy per time.
\section restrictions
keyword 'all' reserved in element_set name
\section related
see \ref man_remove_source
\section default
none
*/
else
if
(
strcmp
(
arg
[
1
],
"source"
)
==
0
)
{
argIdx
=
2
;
get_field
(
arg
,
argIdx
,
thisField
,
thisIndex
);
string
esetName
(
arg
[
argIdx
++
]);
XT_Function
*
f
=
NULL
;
// parse constant
if
(
narg
==
argIdx
+
1
)
{
f
=
XT_Function_Mgr
::
instance
()
->
get_constant_function
(
atof
(
arg
[
argIdx
]));
}
// parse function
else
{
f
=
XT_Function_Mgr
::
instance
()
->
get_function
(
&
(
arg
[
argIdx
]),
narg
-
argIdx
);
}
prescribedDataMgr_
->
fix_source
(
esetName
,
thisField
,
thisIndex
,
f
);
fieldMask_
(
thisField
,
PRESCRIBED_SOURCE
)
=
true
;
match
=
true
;
}
/*! \page man_remove_source fix_modify AtC transfer remove_source
\section syntax
fix_modify AtC transfer remove_source <field> <element_set>
- <field> = field name valid for type of physics
- <element_set> = name of set of elements
\section examples
<TT> fix_modify atc transfer remove_source temperature groupNAME </TT>
\section description
Remove a domain source.
\section restrictions
keyword 'all' reserved in element_set name
\section related
see \ref man_source
\section default
*/
else
if
(
strcmp
(
arg
[
1
],
"remove_source"
)
==
0
)
{
argIdx
=
2
;
get_field
(
arg
,
argIdx
,
thisField
,
thisIndex
);
string
esetName
(
arg
[
argIdx
++
]);
prescribedDataMgr_
->
unfix_source
(
esetName
,
thisField
,
thisIndex
);
fieldMask_
(
thisField
,
PRESCRIBED_SOURCE
)
=
false
;
match
=
true
;
}
/*! \page man_fix_flux fix_modify AtC transfer fix_flux
\section syntax
fix_modify AtC transfer fix_flux <field> <face_set> <value | function>
- <field> = field name valid for type of physics, temperature | electron_temperature
- <face_set> = name of set of element faces
\section examples
<TT> fix_modify atc transfer fix_flux temperature faceSet 10.0 </TT> \n
\section description
Command for fixing normal fluxes e.g. heat_flux.
This command only prescribes the normal component of the physical flux, e.g. heat (energy) flux.
The units are in AtC units, i.e. derived from the LAMMPS length, time, and mass scales.
\section restrictions
Only normal fluxes (Neumann data) can be prescribed.
\section related
see \ref man_unfix_flux
\section default
*/
else
if
(
strcmp
(
arg
[
1
],
"fix_flux"
)
==
0
)
{
argIdx
=
2
;
get_field
(
arg
,
argIdx
,
thisField
,
thisIndex
);
string
fsetName
(
arg
[
argIdx
++
]);
XT_Function
*
f
=
NULL
;
// parse constant
if
(
narg
==
argIdx
+
1
)
{
f
=
XT_Function_Mgr
::
instance
()
->
get_constant_function
(
atof
(
arg
[
argIdx
]));
}
// parse function
else
{
f
=
XT_Function_Mgr
::
instance
()
->
get_function
(
&
(
arg
[
argIdx
]),
narg
-
argIdx
);
}
prescribedDataMgr_
->
fix_flux
(
fsetName
,
thisField
,
thisIndex
,
f
);
fieldMask_
(
thisField
,
PRESCRIBED_SOURCE
)
=
true
;
match
=
true
;
}
/*! \page man_unfix_flux fix_modify AtC transfer unfix_flux
\section syntax
fix_modify AtC transfer fix_flux <field> <face_set> <value | function>
- <field> = field name valid for type of physics, temperature | electron_temperature
- <face_set> = name of set of element faces
\section examples
<TT> fix_modify atc transfer unfix_flux temperature faceSet </TT> \n
\section description
Command for removing prescribed normal fluxes e.g. heat_flux, stress.
\section restrictions
\section related
see \ref man_unfix_flux
\section default
*/
else
if
(
strcmp
(
arg
[
1
],
"unfix_flux"
)
==
0
)
{
argIdx
=
2
;
get_field
(
arg
,
argIdx
,
thisField
,
thisIndex
);
string
fsetName
(
arg
[
argIdx
++
]);
prescribedDataMgr_
->
unfix_flux
(
fsetName
,
thisField
,
thisIndex
);
fieldMask_
(
thisField
,
PRESCRIBED_SOURCE
)
=
false
;
match
=
true
;
}
// pass off to time filter
else
if
(
strcmp
(
arg
[
1
],
"filter"
)
==
0
)
{
match
=
timeFilterManager_
.
modify
(
narg
-
2
,
&
arg
[
2
]);
if
(
match
&&
timeIntegrator_
)
timeIntegrator_
->
force_reset
();
}
}
// pass off to extrinsic model
else
if
(
strcmp
(
arg
[
0
],
"extrinsic"
)
==
0
)
{
match
=
extrinsicModelManager_
.
modify
(
narg
-
1
,
&
arg
[
1
]);
}
/*! \page man_atom_element_map fix_modify AtC transfer atom_element_map
\section syntax
fix_modify AtC transfer atom_element_map <eulerian|lagrangian> <frequency> \n
- frequency (int) : frequency of updating atom-to-continuum maps based on the
current configuration - only for eulerian
\section examples
<TT> fix_modify atc transfer atom_element_map eulerian 100 </TT>
\section description
Changes frame of reference from eulerian to lagrangian and sets the
frequency for which the map from atoms to elements is reformed and
all the attendant data is recalculated.
\section restrictions
Cannot change map type after initialization.
\section related
\section default
lagrangian
*/
else
if
(
strcmp
(
arg
[
0
],
"atom_element_map"
)
==
0
)
{
if
(
strcmp
(
arg
[
1
],
"eulerian"
)
==
0
)
{
atomToElementMapType_
=
EULERIAN
;
atomToElementMapFrequency_
=
atoi
(
arg
[
2
]);
}
else
{
atomToElementMapType_
=
LAGRANGIAN
;
atomToElementMapFrequency_
=
0
;
}
match
=
true
;
}
/*! \page man_neighbor_reset_frequency fix_modify AtC transfer neighbor_reset_frequency
\section syntax
fix_modify AtC transfer neighbor_reset_frequency <frequency>
- frequency (int) : reset neighbor map for AtC analysis every "frequency" steps
\section examples
<TT> fix_modify AtC transfer neighbor_reset_frequency 1000 </TT>
\section description
This command allows the control of how often the lists
of which atoms pairs are in the support of a node
are reformed
\section restrictions
Only currently relevant to Hardy post-processing see \ref man_fix_atc
\section related
This command is different from \ref man_atom_element_map
in that it only concerns neighbor lists internal to ATC not
the map of atom locations to elements which define the shape function
values at atoms.
\section default
The default is for the internal neighbor lists to be recalculated on
the next output step (or sample step if time filtering is used).
*/
else
if
(
strcmp
(
arg
[
0
],
"neighbor_reset_frequency"
)
==
0
)
{
neighborResetFrequency_
=
atoi
(
arg
[
1
]);
match
=
true
;
}
/*! \page man_read_restart fix_modify AtC transfer read_restart
\section syntax
fix_modify AtC transfer read_restart [file_name] \n
\section examples
<TT> fix_modify AtC transfer read_restart ATC_state </TT> \n
\section description
Reads the current state of the fields from a named text-based restart
file.
\section restrictions
The restart file only contains fields and their time derivatives.
The reference positions of the atoms and the commands that initialize
the fix are not saved e.g. an identical mesh containing the same atoms
will have to be recreated.
\section related
see write_restart \ref man_write_restart
\section default
none
*/
else
if
(
strcmp
(
arg
[
0
],
"read_restart"
)
==
0
)
{
restartFileName_
=
arg
[
1
];
useRestart_
=
true
;
match
=
true
;
}
/*! \page man_write_restart fix_modify AtC transfer write_restart
\section syntax
fix_modify AtC transfer write_restart [file_name] \n
\section examples
<TT> fix_modify AtC transfer write_restart restart.mydata </TT> \n
\section description
Dumps the current state of the fields to a named text-based restart file.
This done when the command is invoked and not repeated, unlike the
similar lammps command.
\section restrictions
The restart file only contains fields and their time derivatives.
The reference positions of the atoms and the commands that initialize
the fix are not saved e.g. an identical mesh containing the same atoms
will have to be recreated.
\section related
see read_restart \ref man_read_restart
\section default
none
*/
else
if
(
strcmp
(
arg
[
0
],
"write_restart"
)
==
0
)
{
string
restartFileName
(
arg
[
1
]);
OUTPUT_LIST
data
;
write_restart_data
(
restartFileName
,
data
);
match
=
true
;
}
// pass off to fe engine
if
(
!
match
)
{
match
=
feEngine_
->
modify
(
narg
,
arg
);
}
return
match
;
}
//--------------------------------------------------
// parse_boundary_integration
// parses the boundary integration to determine
// the type of boundary integration being used
//--------------------------------------------------
ATC
::
ATC_Transfer
::
BoundaryIntegrationType
ATC_Transfer
::
parse_boundary_integration
(
int
narg
,
char
**
arg
,
const
set
<
pair
<
int
,
int
>
>
*
boundaryFaceSet
)
{
int
argIndex
=
0
;
BoundaryIntegrationType
myBoundaryIntegrationType
;
if
(
narg
>
0
)
{
if
(
strcmp
(
arg
[
argIndex
],
"faceset"
)
==
0
)
{
argIndex
++
;
myBoundaryIntegrationType
=
FE_QUADRATURE
;
string
name
(
arg
[
argIndex
]);
boundaryFaceSet
=
&
(
(
feEngine_
->
get_feMesh
())
->
get_faceset
(
name
));
set_boundary_face_set
(
boundaryFaceSet
);
}
else
if
(
strcmp
(
arg
[
argIndex
],
"interpolate"
)
==
0
)
{
myBoundaryIntegrationType
=
FE_INTERPOLATION
;
}
else
if
(
strcmp
(
arg
[
argIndex
],
"no_boundary"
)
==
0
)
{
myBoundaryIntegrationType
=
NO_QUADRATURE
;
}
else
{
throw
ATC_Error
(
0
,
"Bad boundary integration type"
);
}
}
else
{
// default is interpolation
myBoundaryIntegrationType
=
FE_INTERPOLATION
;
}
set_boundary_integration_type
(
myBoundaryIntegrationType
);
return
myBoundaryIntegrationType
;
}
//--------------------------------------------------
// helper function for parser
// handles : "displacement x" and "temperature" by indexing argIdx
// for fluxes : only normal fluxes can be prescribed
//--------------------------------------------------
void
ATC_Transfer
::
get_field
(
/*const*/
char
**
args
,
int
&
argIdx
,
FieldName
&
thisField
,
int
&
thisIndex
)
{
string
thisName
=
args
[
argIdx
++
];
thisField
=
string_to_field
(
thisName
);
map
<
FieldName
,
int
>::
const_iterator
iter
=
fieldSizes_
.
find
(
thisField
);
if
(
iter
==
fieldSizes_
.
end
())
{
throw
ATC_Error
(
0
,
"Bad field name: "
+
thisName
);
}
string
thisDim
=
args
[
argIdx
];
thisIndex
=
0
;
if
(
string_to_index
(
thisDim
,
thisIndex
))
{
if
(
!
(
thisIndex
<
fieldSizes_
[
thisField
])
)
{
throw
ATC_Error
(
0
,
"Bad field dimension "
+
thisDim
);
}
argIdx
++
;
}
}
void
ATC_Transfer
::
init_filter
()
{
if
(
timeFilterManager_
.
need_reset
())
{
map
<
FieldName
,
int
>::
const_iterator
field
;
for
(
field
=
fieldSizes_
.
begin
();
field
!=
fieldSizes_
.
end
();
field
++
)
{
FieldName
thisField
=
field
->
first
;
int
thisSize
=
field
->
second
;
// Allocate restricted fields
fieldRateNdFiltered_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
fieldRateNdOld_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
dot_fieldRateNdFiltered_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
dot_fieldRateNdOld_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
}
}
}
void
ATC_Transfer
::
set_xref
()
{
double
**
x
=
lammpsInterface_
->
xatom
();
for
(
int
i
=
0
;
i
<
lammpsInterface_
->
nmax
();
i
++
)
{
for
(
int
j
=
0
;
j
<
nsd_
;
j
++
)
{
xref_
[
i
][
j
]
=
x
[
i
][
j
];
}
}
Xprd_
=
lammpsInterface_
->
domain_xprd
();
Yprd_
=
lammpsInterface_
->
domain_yprd
();
Zprd_
=
lammpsInterface_
->
domain_zprd
();
XY_
=
lammpsInterface_
->
domain_xy
();
XZ_
=
lammpsInterface_
->
domain_xz
();
YZ_
=
lammpsInterface_
->
domain_yz
();
}
void
ATC_Transfer
::
initialize
()
{
// initialized_ is set to true by derived class initialize()
int
i
,
j
;
if
(
!
initialized_
)
{
if
(
!
trackCharge_
)
{
fieldSizes_
.
erase
(
CHARGE_DENSITY
);
prescribedDataMgr_
->
remove_field
(
CHARGE_DENSITY
);
for
(
i
=
0
;
i
<
NUM_FLUX
;
i
++
)
fieldMask_
(
CHARGE_DENSITY
,
i
)
=
false
;
}
// Initialize finite elements shape function, mapping variables
feEngine_
->
initialize
();
nNodes_
=
feEngine_
->
get_nNodes
();
nsd_
=
feEngine_
->
get_nsd
();
if
(
nsd_
!=
lammpsInterface_
->
dimension
())
throw
ATC_Error
(
0
,
"Spatial dimensions inconsistent between LAMMPS and ATC"
);
// set reference position of atoms
set_xref
();
}
// determine which elements are strictly in the MD region and mask out for intrinsic fields
// start by finding local atom counts
reset_nlocal
();
int
nelts
=
feEngine_
->
get_feMesh
()
->
get_nElements
();
elementMask_
.
reset
(
nelts
);
if
(
atomQuadForInternal_
)
{
elementMask_
=
true
;
}
else
{
// determine which elements contain internal atoms
Array
<
int
>
hasInternal
(
nelts
);
Array
<
int
>
hasInternal_local
(
nelts
);
hasInternal_local
=
0
;
for
(
int
i
=
0
;
i
<
nLocal_
;
++
i
)
{
DENS_VEC
coords
(
3
);
coords
(
0
)
=
atomicCoords_
(
0
,
i
);
coords
(
1
)
=
atomicCoords_
(
1
,
i
);
coords
(
2
)
=
atomicCoords_
(
2
,
i
);
int
eltID
=
feEngine_
->
get_feMesh
()
->
map_to_element
(
coords
);
hasInternal_local
(
eltID
)
=
1
;
}
// swap contributions
lammpsInterface_
->
logical_or
(
hasInternal_local
.
get_ptr
(),
hasInternal
.
get_ptr
(),
hasInternal
.
size
());
// determine which elements contain ghost atoms
Array
<
int
>
hasGhost
(
nelts
);
Array
<
int
>
hasGhost_local
(
nelts
);
hasGhost_local
=
0
;
for
(
int
i
=
0
;
i
<
nLocalGhost_
;
++
i
)
{
DENS_VEC
coords
(
3
);
coords
(
0
)
=
ghostAtomCoords_
(
0
,
i
);
coords
(
1
)
=
ghostAtomCoords_
(
1
,
i
);
coords
(
2
)
=
ghostAtomCoords_
(
2
,
i
);
int
eltID
=
feEngine_
->
get_feMesh
()
->
map_to_element
(
coords
);
hasGhost_local
(
eltID
)
=
1
;
}
//swap contributions
lammpsInterface_
->
logical_or
(
hasGhost_local
.
get_ptr
(),
hasGhost
.
get_ptr
(),
hasGhost
.
size
());
for
(
int
i
=
0
;
i
<
nelts
;
++
i
)
elementMask_
(
i
)
=
!
hasInternal
(
i
)
||
hasGhost
(
i
);
}
set
<
int
>
&
nullElements
=
feEngine_
->
null_elements
();
set
<
int
>::
const_iterator
iset
;
for
(
iset
=
nullElements
.
begin
();
iset
!=
nullElements
.
end
();
iset
++
)
{
int
ielem
=
*
iset
;
elementMask_
(
ielem
)
=
false
;
}
// set atomic weights
reset_shape_functions
();
// creates atom->element map
set_atomic_weights
();
// group atoms by material
int
nMatls
=
physicsModel_
->
get_nMaterials
();
atomMaterialGroups_
.
reset
(
nMatls
);
for
(
int
i
=
0
;
i
<
nLocal_
;
i
++
)
{
int
matId
=
elementToMaterialMap_
(
atomToElementMap_
(
i
));
atomMaterialGroups_
(
matId
).
insert
(
i
);
// not internalToAtom_
}
if
(
atomToElementMapType_
==
EULERIAN
&&
nMatls
>
1
)
{
throw
ATC_Error
(
0
,
" multiple materials & eulerian map not supported"
);
}
// compute necessary restructuring matrices
maskMat_
.
reset
(
nNodes_
,
nNodes_
);
if
(
!
useLocalizedLambda_
)
{
for
(
i
=
0
;
i
<
nNodes_
;
++
i
)
if
(
shpWeight_
(
i
,
i
)
>
0.
)
maskMat_
(
i
,
i
)
=
1.
;
}
else
{
// select nodes for thermostat
// 1) their shape functions have both mask and ghost atoms in their support
DENS_VEC
hasInternalLocal
(
nNodes_
);
DENS_VEC
hasInternal
(
nNodes_
);
if
(
nLocal_
>
0
)
hasInternalLocal
=
shpFcn_
.
col_sum
();
lammpsInterface_
->
allsum
(
hasInternalLocal
.
get_ptr
(),
hasInternal
.
get_ptr
(),
nNodes_
);
DENS_VEC
hasGhostLocal
(
nNodes_
);
DENS_VEC
hasGhost
(
nNodes_
);
if
(
nLocalGhost_
>
0
)
hasGhostLocal
=
shpFcnGhost_
.
col_sum
();
lammpsInterface_
->
allsum
(
hasGhostLocal
.
get_ptr
(),
hasGhost
.
get_ptr
(),
nNodes_
);
for
(
i
=
0
;
i
<
nNodes_
;
++
i
)
{
if
(
hasInternal
(
i
)
>
0.
&&
hasGhost
(
i
)
>
0.
)
{
maskMat_
(
i
,
i
)
=
1.
;
}
}
// 2) they are in a specified boundary faceset
set
<
string
>::
const_iterator
iter
;
for
(
iter
=
boundaryFaceNames_
.
begin
();
iter
!=
boundaryFaceNames_
.
end
();
iter
++
)
{
set
<
int
>
nodeSet
;
feEngine_
->
get_feMesh
()
->
faceset_to_nodeset
(
*
iter
,
nodeSet
);
set
<
int
>::
const_iterator
nodeIter
;
for
(
nodeIter
=
nodeSet
.
begin
();
nodeIter
!=
nodeSet
.
end
();
nodeIter
++
)
{
maskMat_
(
*
nodeIter
,
*
nodeIter
)
=
1.
;
}
}
}
// determine non-zero values and row-column indices
// also used in creating overlap map and mapping from atoms to elements
// note we assume the overlap map is time-independent
nodeToOverlapMap_
.
reset
(
nNodes_
);
nNodeOverlap_
=
0
;
for
(
i
=
0
;
i
<
nNodes_
;
i
++
)
{
if
(
maskMat_
(
i
,
i
)
>
0.
)
{
nodeToOverlapMap_
(
i
)
=
nNodeOverlap_
;
nNodeOverlap_
++
;
}
else
{
nodeToOverlapMap_
(
i
)
=
-
1
;
}
}
overlapToNodeMap_
.
reset
(
nNodeOverlap_
);
int
counter
=
0
;
for
(
i
=
0
;
i
<
nNodes_
;
i
++
)
{
if
(
nodeToOverlapMap_
(
i
)
!=-
1
)
{
overlapToNodeMap_
(
counter
)
=
i
;
counter
++
;
}
}
// reset local overlap shape function storage
reset_NhatOverlap
();
// Shape sparse matrix used in Hoover-type thermostats
// first get local pattern sized nNodeOverlap X nNodeOverlap
//compute_consistent_md_mass_matrix(NhatOverlap_,tmp);
DENS_MAT
tmpLocal
(
nNodeOverlap_
,
nNodeOverlap_
);
DENS_MAT
tmp
(
nNodeOverlap_
,
nNodeOverlap_
);
if
(
nLocalLambda_
>
0
)
{
DENS_VEC
T
(
nLocalLambda_
);
T
=
1.0
;
for
(
i
=
0
;
i
<
nNodeOverlap_
;
i
++
)
{
for
(
j
=
0
;
j
<
nNodeOverlap_
;
j
++
)
{
for
(
int
iatom
=
0
;
iatom
<
nLocalLambda_
;
iatom
++
)
{
tmpLocal
(
i
,
j
)
+=
NhatOverlap_
(
iatom
,
i
)
*
T
(
iatom
)
*
NhatOverlap_
(
iatom
,
j
);
}
}
}
}
// second accumulate total pattern
lammpsInterface_
->
allsum
(
tmpLocal
.
get_ptr
(),
tmp
.
get_ptr
(),
tmp
.
size
());
// third extract non-zero entries & construct
M_T_Template
.
reset
(
nNodeOverlap_
,
nNodeOverlap_
);
for
(
i
=
0
;
i
<
nNodeOverlap_
;
i
++
)
{
for
(
j
=
0
;
j
<
nNodeOverlap_
;
j
++
)
{
if
(
abs
(
tmp
(
i
,
j
))
>
0
)
{
M_T_Template
.
add
(
i
,
j
,
0.
);
}
}
}
M_T_Template
.
compress
();
// \int_\Omega N_I dV
NodeVolumes_
.
reset
(
nNodes_
,
nNodes_
);
invNodeVolumes_
.
reset
(
nNodes_
,
nNodes_
);
feEngine_
->
compute_lumped_mass_matrix
(
NodeVolumes_
);
invNodeVolumes_
=
NodeVolumes_
.
inv
();
// this is exact only for uniform meshes and certain types of atomic weights
// \int_{\Omega_MD} N_I dV = \sum_\alpha N_I\alpha V_\alpha
// DENS_VEC atom_weights_vec(atomicWeights_.size(),
// atomicWeights_.get_ptr());
DENS_VEC
subsetNodeVolume
(
NodeVolumes_
.
size
());
DENS_VEC
unitVec
(
atomicWeights_
.
size
());
unitVec
=
1.0
;
restrict_unscaled
(
unitVec
,
subsetNodeVolume
);
fluxMask_
.
reset
(
invNodeVolumes_
*
subsetNodeVolume
);
DIAG_MAT
id
(
fluxMask_
.
nRows
(),
fluxMask_
.
nCols
());
id
=
1.0
;
fluxMaskComplement_
=
id
+
-
1.0
*
fluxMask_
;
// set flux masks for nodes we can tell by geometry
DENS_VEC
hasInternalLocal
(
nNodes_
);
DENS_VEC
hasInternal
(
nNodes_
);
if
(
nLocal_
>
0
)
hasInternalLocal
=
shpFcn_
.
col_sum
();
lammpsInterface_
->
allsum
(
hasInternalLocal
.
get_ptr
(),
hasInternal
.
get_ptr
(),
nNodes_
);
DENS_VEC
hasGhostLocal
(
nNodes_
);
DENS_VEC
hasGhost
(
nNodes_
);
if
(
nLocalGhost_
>
0
)
hasGhostLocal
=
shpFcnGhost_
.
col_sum
();
lammpsInterface_
->
allsum
(
hasGhostLocal
.
get_ptr
(),
hasGhost
.
get_ptr
(),
nNodes_
);
for
(
i
=
0
;
i
<
nNodes_
;
++
i
)
{
if
(
hasInternal
(
i
)
>
0.
)
{
if
(
hasGhost
(
i
)
==
0.
)
{
fluxMask_
(
i
,
i
)
=
1.
;
fluxMaskComplement_
(
i
,
i
)
=
0.
;
}
}
else
{
fluxMask_
(
i
,
i
)
=
0.
;
fluxMaskComplement_
(
i
,
i
)
=
1.
;
}
}
// Initialize fields and set initial conditions
if
(
!
initialized_
)
{
nodeType_
.
reset
(
nNodes_
);
for
(
int
iNode
=
0
;
iNode
<
nNodes_
;
++
iNode
)
nodeType_
(
iNode
)
=
check_shape_function_type
(
iNode
);
map
<
FieldName
,
int
>::
const_iterator
field
;
for
(
field
=
fieldSizes_
.
begin
();
field
!=
fieldSizes_
.
end
();
field
++
)
{
FieldName
thisField
=
field
->
first
;
int
thisSize
=
field
->
second
;
// Allocate fields, initialize to default values, set up initial schedule
fields_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
dot_fields_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
ddot_fields_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
dddot_fields_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
dot_fieldsMD_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
ddot_fieldsMD_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
dot_dot_fieldsMD_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
// Allocate restricted fields
fieldRateNdFiltered_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
fieldRateNdOld_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
dot_fieldRateNdFiltered_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
dot_fieldRateNdOld_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
fieldNdFiltered_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
fieldNdOld_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
// Allocate auxilliary storage
auxStorage_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
// Dimension finite element rhs matrix
rhs_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
rhsAtomDomain_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
sources_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
extrinsicSources_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
boundaryFlux_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
if
(
is_intrinsic
(
thisField
))
{
// compute FE mass matrix in full domain
massMatInv_
[
thisField
].
reset
(
nNodes_
,
nNodes_
);
Array
<
FieldName
>
massMask
(
1
);
massMask
(
0
)
=
thisField
;
feEngine_
->
compute_lumped_mass_matrix
(
massMask
,
fields_
,
physicsModel_
,
elementToMaterialMap_
,
massMats_
,
&
elementMask_
);
// atomic quadrature for FE mass matrix in atomic domain
map
<
FieldName
,
DIAG_MAT
>
massMatAq
;
feEngine_
->
compute_lumped_mass_matrix
(
massMask
,
fields_
,
physicsModel_
,
atomMaterialGroups_
,
atomicWeightsMask_
,
shpFcnMasked_
,
massMatAq
);
// remove contributions from overlap by approximate quadrature
// and fully remove contributions from internal
massMats_
[
thisField
]
-=
massMatAq
[
thisField
];
for
(
int
iNode
=
0
;
iNode
<
nNodes_
;
iNode
++
)
if
(
nodeType_
(
iNode
)
==
MD_ONLY
)
massMats_
[
thisField
](
iNode
,
iNode
)
=
0.
;
// set up mass MD matrices
massMatMDInv_
[
thisField
].
reset
(
nNodes_
,
nNodes_
);
compute_md_mass_matrix
(
thisField
,
massMatsMD_
);
massMats_
[
thisField
]
+=
massMatsMD_
[
thisField
];
// compute inverse mass matrices since we're using lumped masses
for
(
int
iNode
=
0
;
iNode
<
nNodes_
;
iNode
++
)
{
if
(
fabs
(
massMatsMD_
[
thisField
](
iNode
,
iNode
))
>
0
)
{
massMatMDInv_
[
thisField
](
iNode
,
iNode
)
=
1.
/
massMatsMD_
[
thisField
](
iNode
,
iNode
);
}
if
(
fabs
(
massMats_
[
thisField
](
iNode
,
iNode
))
>
0
)
{
massMatInv_
[
thisField
](
iNode
,
iNode
)
=
1.
/
massMats_
[
thisField
](
iNode
,
iNode
);
}
}
}
else
{
// no MD mass matrices needed, regular matrices computed in extrinsic model
massMats_
[
thisField
].
reset
(
nNodes_
,
nNodes_
);
massMatInv_
[
thisField
].
reset
(
nNodes_
,
nNodes_
);
}
}
// compute inverse consistent mass for the continuity equation
SPAR_MAT
consistentMass
;
feEngine_
->
compute_mass_matrix
(
consistentMass
);
consistentMass
.
compress
();
consistentMassInverse_
=
inv
(
consistentMass
.
dense_copy
());
// Apply integration masking and new ICs
// initialize schedule derivatives
try
{
set_initial_conditions
();
}
catch
(
ATC
::
ATC_Error
&
atcError
)
{
if
(
!
useRestart_
)
throw
;
}
}
// set lists of fixed nodes per field
map
<
FieldName
,
int
>::
const_iterator
field
;
for
(
field
=
fieldSizes_
.
begin
();
field
!=
fieldSizes_
.
end
();
field
++
)
{
FieldName
thisField
=
field
->
first
;
int
thisSize
=
field
->
second
;
isFixedNode_
[
thisField
].
reset
(
nNodes_
,
thisSize
);
for
(
int
i
=
0
;
i
<
nNodes_
;
i
++
)
for
(
int
j
=
0
;
j
<
thisSize
;
j
++
)
isFixedNode_
[
thisField
](
i
,
j
)
=
prescribedDataMgr_
->
is_fixed
(
i
,
thisField
,
j
);
}
return
;
}
void
ATC_Transfer
::
pre_final_integrate
()
{
// atomSwtich is triggered by a lammps call to pre_exchange
// (see ATC_Transfer.h)
int
localSwitch
=
0
;
if
(
atomSwitch_
||
resetSwitch_
)
localSwitch
=
1
;
lammpsInterface_
->
int_allmax
(
&
localSwitch
,
&
globalSwitch_
);
if
(
globalSwitch_
>
0
)
{
reset_nlocal
();
reset_shape_functions
();
set_atomic_weights
();
reset_NhatOverlap
();
atomSwitch_
=
false
;
globalSwitch_
=
0
;
}
else
if
(
atomToElementMapType_
==
EULERIAN
&&
stepCounter_
%
atomToElementMapFrequency_
==
0
)
{
reset_coordinates
();
reset_shape_functions
();
set_atomic_weights
();
reset_NhatOverlap
();
}
}
void
ATC_Transfer
::
init_integrate_velocity
()
{
int
nlocal
=
lammpsInterface_
->
nlocal
();
int
*
type
=
lammpsInterface_
->
atom_type
();
int
*
mask
=
lammpsInterface_
->
atom_mask
();
double
*
mass
=
lammpsInterface_
->
atom_mass
();
double
**
v
=
lammpsInterface_
->
vatom
();
double
**
f
=
lammpsInterface_
->
fatom
();
double
dtf
=
0.5
*
lammpsInterface_
->
dt
()
*
lammpsInterface_
->
ftm2v
();
if
(
mass
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit_
)
{
double
dtfm
=
dtf
/
mass
[
type
[
i
]];
v
[
i
][
0
]
+=
dtfm
*
f
[
i
][
0
];
v
[
i
][
1
]
+=
dtfm
*
f
[
i
][
1
];
v
[
i
][
2
]
+=
dtfm
*
f
[
i
][
2
];
}
}
}
else
{
double
*
rmass
=
lammpsInterface_
->
atom_rmass
();
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit_
)
{
double
dtfm
=
dtf
/
rmass
[
i
];
v
[
i
][
0
]
+=
dtfm
*
f
[
i
][
0
];
v
[
i
][
1
]
+=
dtfm
*
f
[
i
][
1
];
v
[
i
][
2
]
+=
dtfm
*
f
[
i
][
2
];
}
}
}
}
void
ATC_Transfer
::
init_integrate_position
()
{
int
nlocal
=
lammpsInterface_
->
nlocal
();
double
**
x
=
lammpsInterface_
->
xatom
();
double
**
v
=
lammpsInterface_
->
vatom
();
double
dtv
=
lammpsInterface_
->
dt
();
int
*
mask
=
lammpsInterface_
->
atom_mask
();
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit_
)
{
x
[
i
][
0
]
+=
dtv
*
v
[
i
][
0
];
x
[
i
][
1
]
+=
dtv
*
v
[
i
][
1
];
x
[
i
][
2
]
+=
dtv
*
v
[
i
][
2
];
}
}
}
void
ATC_Transfer
::
final_integrate
()
{
int
nlocal
=
lammpsInterface_
->
nlocal
();
int
*
type
=
lammpsInterface_
->
atom_type
();
int
*
mask
=
lammpsInterface_
->
atom_mask
();
double
*
mass
=
lammpsInterface_
->
atom_mass
();
double
**
v
=
lammpsInterface_
->
vatom
();
double
**
f
=
lammpsInterface_
->
fatom
();
double
dtf
=
0.5
*
lammpsInterface_
->
dt
()
*
lammpsInterface_
->
ftm2v
();
if
(
mass
)
{
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit_
)
{
double
dtfm
=
dtf
/
mass
[
type
[
i
]];
v
[
i
][
0
]
+=
dtfm
*
f
[
i
][
0
];
v
[
i
][
1
]
+=
dtfm
*
f
[
i
][
1
];
v
[
i
][
2
]
+=
dtfm
*
f
[
i
][
2
];
}
}
}
else
{
double
*
rmass
=
lammpsInterface_
->
atom_rmass
();
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
{
if
(
mask
[
i
]
&
groupbit_
)
{
double
dtfm
=
dtf
/
rmass
[
i
];
v
[
i
][
0
]
+=
dtfm
*
f
[
i
][
0
];
v
[
i
][
1
]
+=
dtfm
*
f
[
i
][
1
];
v
[
i
][
2
]
+=
dtfm
*
f
[
i
][
2
];
}
}
}
}
void
ATC_Transfer
::
finish
()
{
// possibly resetting flags and deallocating corresponding memory
// FE Engine
feEngine_
->
finish
();
return
;
}
//--------------------------------------------------------------
// create_physics_model
// - method to create physics model
//--------------------------------------------------------------
void
ATC_Transfer
::
create_physics_model
(
const
PhysicsType
&
physicsType
,
string
matFileName
)
{
if
(
physicsModel_
)
{
throw
ATC_Error
(
0
,
"Attempted to create PhysicsModel multiple times in ATC_Transfer"
);
}
// Create PhysicsModel based on physicsType
switch
(
physicsType
)
{
case
NO_PHYSICS
:
break
;
case
THERMAL
:
physicsModel_
=
new
PhysicsModelThermal
(
matFileName
,
this
);
break
;
default
:
throw
ATC_Error
(
0
,
"Unknown physics type in ATC_Transfer::create_physics_model"
);
}
}
//--------------------------------------------------------------
/** method to trigger construction of mesh data after mesh construction */
//--------------------------------------------------------------
void
ATC_Transfer
::
initialize_mesh_data
(
void
)
{
int
nelts
=
feEngine_
->
get_feMesh
()
->
get_nElements
();
elementToMaterialMap_
.
reset
(
nelts
);
elementToMaterialMap_
=
0
;
construct_prescribed_data_manager
();
}
//--------------------------------------------------------------
/** method to add new fields to the included list */
//--------------------------------------------------------------
void
ATC_Transfer
::
add_fields
(
map
<
FieldName
,
int
>
&
newFieldSizes
)
{
map
<
FieldName
,
int
>::
const_iterator
field
;
for
(
field
=
newFieldSizes
.
begin
();
field
!=
newFieldSizes
.
end
();
field
++
)
{
FieldName
thisField
=
field
->
first
;
int
thisSize
=
field
->
second
;
if
(
fieldSizes_
.
find
(
thisField
)
==
fieldSizes_
.
end
())
{
fieldSizes_
[
thisField
]
=
thisSize
;
}
}
}
};
// namespace ATC
Event Timeline
Log In to Comment