Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73381551
FE_Engine.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
Sun, Jul 21, 00:03
Size
108 KB
Mime Type
text/x-c
Expires
Tue, Jul 23, 00:03 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
19194973
Attached To
rLAMMPS lammps
FE_Engine.cpp
View Options
#include "FE_Engine.h"
#include "ATC_Transfer.h"
#include "FE_Element.h"
#include "Function.h"
#include "PhysicsModel.h"
#include "KernelFunction.h"
#include "Utility.h"
#include "MPI_Wrappers.h"
#include <stdio.h>
#include <stdlib.h>
#include <map>
#include <string>
using
namespace
std
;
using
ATC_Utility
::
is_numeric
;
using
ATC_Utility
::
to_string
;;
using
MPI_Wrappers
::
allsum
;
using
MPI_Wrappers
::
int_allsum
;
using
MPI_Wrappers
::
rank_zero
;
using
MPI_Wrappers
::
print_msg
;
using
MPI_Wrappers
::
print_msg_once
;
namespace
ATC
{
static
const
double
tol_sparse
=
1.e-30
;
//tolerance for compaction from dense
//-----------------------------------------------------------------
FE_Engine
::
FE_Engine
(
MPI_Comm
comm
)
:
communicator_
(
comm
),
feMesh_
(
NULL
),
initialized_
(
false
),
outputManager_
()
{
// Nothing to do here
}
//-----------------------------------------------------------------
FE_Engine
::~
FE_Engine
()
{
if
(
feMesh_
)
delete
feMesh_
;
}
//-----------------------------------------------------------------
void
FE_Engine
::
initialize
()
{
if
(
!
feMesh_
)
throw
ATC_Error
(
"FE_Engine has no mesh"
);
if
(
!
initialized_
)
{
// set up work spaces
nNodesPerElement_
=
feMesh_
->
num_nodes_per_element
();
nIPsPerElement_
=
feMesh_
->
num_ips_per_element
();
nIPsPerFace_
=
feMesh_
->
num_ips_per_face
();
nSD_
=
feMesh_
->
num_spatial_dimensions
();
nElems_
=
feMesh_
->
num_elements
();
nNodesUnique_
=
feMesh_
->
num_nodes_unique
();
nNodes_
=
feMesh_
->
num_nodes
();
// arrays & matrices
_weights_
.
reset
(
nIPsPerElement_
,
nIPsPerElement_
);
_N_
.
reset
(
nIPsPerElement_
,
nNodesPerElement_
);
_dN_
.
assign
(
nSD_
,
DENS_MAT
(
nIPsPerElement_
,
nNodesPerElement_
));
_Nw_
.
reset
(
nIPsPerElement_
,
nNodesPerElement_
);
_dNw_
.
assign
(
nSD_
,
DENS_MAT
(
nIPsPerElement_
,
nNodesPerElement_
));
_Bfluxes_
.
assign
(
nSD_
,
DENS_MAT
());
// faces
_fweights_
.
reset
(
nIPsPerElement_
,
nIPsPerElement_
);
_fN_
.
reset
(
nIPsPerFace_
,
nNodesPerElement_
);
_fdN_
.
assign
(
nSD_
,
DENS_MAT
(
nIPsPerFace_
,
nNodesPerElement_
));
_nN_
.
assign
(
nSD_
,
DENS_MAT
(
nIPsPerFace_
,
nNodesPerElement_
));
// remove specified elements
if
(
nullElements_
.
size
()
>
0
)
delete_elements
(
nullElements_
);
initialized_
=
true
;
}
}
void
FE_Engine
::
partition_mesh
()
{
if
(
is_partitioned
())
return
;
feMesh_
->
partition_mesh
();
// now do all FE_Engine data structure partitioning
// partition nullElements_
/*for (vector<int>::iterator elemsIter = feMesh_->processor_elts().begin();
elemsIter != feMesh_->processor_elts().end();
++elemsIter)
{
if (nullElements_.find(*elemsIter) != nullElements_.end()) {
myNullElements_.insert(map_elem_to_myElem(*elemsIter));
}
}*/
}
void
FE_Engine
::
departition_mesh
()
{
if
(
!
is_partitioned
())
return
;
feMesh_
->
departition_mesh
();
}
void
FE_Engine
::
set_quadrature
(
FeIntQuadrature
type
,
bool
temp
)
const
{
if
(
!
feMesh_
)
throw
ATC_Error
(
"FE_Engine has no mesh"
);
feMesh_
->
set_quadrature
(
type
);
if
(
!
temp
)
quadrature_
=
type
;
int
nIPsPerElement_new
=
feMesh_
->
num_ips_per_element
();
int
nIPsPerFace_new
=
feMesh_
->
num_ips_per_face
();
if
(
nIPsPerElement_
!=
nIPsPerElement_new
)
{
// arrays & matrices
nIPsPerElement_
=
nIPsPerElement_new
;
_weights_
.
resize
(
nIPsPerElement_
,
nIPsPerElement_
);
_N_
.
resize
(
nIPsPerElement_
,
nNodesPerElement_
);
_dN_
.
assign
(
nSD_
,
DENS_MAT
(
nIPsPerElement_
,
nNodesPerElement_
));
_Nw_
.
reset
(
nIPsPerElement_
,
nNodesPerElement_
);
_dNw_
.
assign
(
nSD_
,
DENS_MAT
(
nIPsPerElement_
,
nNodesPerElement_
));
}
if
(
nIPsPerFace_
!=
nIPsPerFace_new
)
{
// faces
nIPsPerFace_
=
nIPsPerFace_new
;
_fweights_
.
reset
(
nIPsPerElement_
,
nIPsPerElement_
);
_fN_
.
reset
(
nIPsPerFace_
,
nNodesPerElement_
);
_fdN_
.
assign
(
nSD_
,
DENS_MAT
(
nIPsPerFace_
,
nNodesPerElement_
));
_nN_
.
assign
(
nSD_
,
DENS_MAT
(
nIPsPerFace_
,
nNodesPerElement_
));
}
}
//-----------------------------------------------------------------
bool
FE_Engine
::
modify
(
int
narg
,
char
**
arg
)
{
bool
match
=
false
;
/*! \page man_mesh_create fix_modify AtC mesh create
\section syntax
fix_modify AtC mesh create <nx> <ny> <nz> <region-id>
<f|p> <f|p> <f|p> \n
- nx ny nz = number of elements in x, y, z
- region-id = id of region that is to be meshed
- f p p = periodicity flags for x, y, z
\section examples
<TT> fix_modify AtC mesh create 10 1 1 feRegion p p p </TT> \n
\section description
Creates a uniform mesh in a rectangular region
\section restrictions
Creates only uniform rectangular grids in a rectangular region
\section related
\ref man_mesh_quadrature
\section default
When created, mesh defaults to gauss2 (2-point Gaussian) quadrature.
Use "mesh quadrature" command to change quadrature style.
*/
int
argIdx
=
0
;
if
(
strcmp
(
arg
[
argIdx
],
"mesh"
)
==
0
)
{
argIdx
++
;
// create mesh
if
(
strcmp
(
arg
[
argIdx
],
"create"
)
==
0
)
{
if
(
feMesh_
)
throw
ATC_Error
(
"FE Engine already has a mesh"
);
argIdx
++
;
int
nx
=
atoi
(
arg
[
argIdx
++
]);
int
ny
=
atoi
(
arg
[
argIdx
++
]);
int
nz
=
atoi
(
arg
[
argIdx
++
]);
string
box
=
arg
[
argIdx
++
];
Array
<
bool
>
periodicity
(
3
);
periodicity
(
0
)
=
(
strcmp
(
arg
[
argIdx
++
],
"p"
)
==
0
)
?
true
:
false
;
periodicity
(
1
)
=
(
strcmp
(
arg
[
argIdx
++
],
"p"
)
==
0
)
?
true
:
false
;
periodicity
(
2
)
=
(
strcmp
(
arg
[
argIdx
++
],
"p"
)
==
0
)
?
true
:
false
;
if
(
argIdx
<
narg
)
{
Array
<
double
>
dx
(
nx
),
dy
(
ny
),
dz
(
nz
);
dx
=
0
;
dy
=
0
;
dz
=
0
;
while
(
argIdx
<
narg
)
{
if
(
strcmp
(
arg
[
argIdx
],
"dx"
)
==
0
)
{
parse_partitions
(
++
argIdx
,
narg
,
arg
,
0
,
dx
);
}
else
if
(
strcmp
(
arg
[
argIdx
],
"dy"
)
==
0
)
{
parse_partitions
(
++
argIdx
,
narg
,
arg
,
1
,
dy
);
}
else
if
(
strcmp
(
arg
[
argIdx
],
"dz"
)
==
0
)
{
parse_partitions
(
++
argIdx
,
narg
,
arg
,
2
,
dz
);
}
}
create_mesh
(
dx
,
dy
,
dz
,
box
.
c_str
(),
periodicity
);
}
else
{
create_mesh
(
nx
,
ny
,
nz
,
box
.
c_str
(),
periodicity
);
}
quadrature_
=
GAUSS2
;
match
=
true
;
}
/*! \page man_mesh_quadrature fix_modify AtC mesh quadrature
\section syntax
fix_modify AtC mesh quadrature <quad>
- quad = one of <nodal|gauss1|gauss2|gauss3|face> --- when a mesh is created it defaults to gauss2, use this call to change it after the fact
\section examples
<TT> fix_modify AtC mesh quadrature face </TT>
\section description
(Re-)assigns the quadrature style for the existing mesh.
\section restrictions
\section related
\ref man_mesh_create
\section default
none
*/
else
if
(
strcmp
(
arg
[
argIdx
],
"quadrature"
)
==
0
)
{
argIdx
++
;
string
quadStr
=
arg
[
argIdx
];
FeIntQuadrature
quadEnum
=
string_to_FIQ
(
quadStr
);
set_quadrature
(
quadEnum
);
match
=
true
;
}
/*! \page man_mesh_read fix_modify AtC mesh read
\section syntax
fix_modify AtC mesh read <filename> <f|p> <f|p> <f|p>
- filename = name of file containing mesh to be read
- f p p = periodicity flags for x, y, z
\section examples
<TT> fix_modify AtC mesh read myComponent.mesh p p p </TT> \n
<TT> fix_modify AtC mesh read myOtherComponent.exo </TT>
\section description
Reads a mesh from a text or exodus file, and assigns periodic
boundary conditions if needed.
\section restrictions
\section related
\section default
periodicity flags are false by default
*/
else
if
(
strcmp
(
arg
[
argIdx
],
"read"
)
==
0
)
{
argIdx
++
;
string
meshFile
=
arg
[
argIdx
++
];
Array
<
bool
>
periodicity
(
3
);
periodicity
=
false
;
if
(
argIdx
<
narg
)
{
periodicity
(
0
)
=
(
strcmp
(
arg
[
argIdx
++
],
"p"
)
==
0
)
?
true
:
false
;
periodicity
(
1
)
=
(
strcmp
(
arg
[
argIdx
++
],
"p"
)
==
0
)
?
true
:
false
;
periodicity
(
2
)
=
(
strcmp
(
arg
[
argIdx
++
],
"p"
)
==
0
)
?
true
:
false
;
}
read_mesh
(
meshFile
,
periodicity
);
if
(
periodicity
(
0
)
||
periodicity
(
1
)
||
periodicity
(
2
))
{
meshFile
=
"periodic_"
+
meshFile
;
stringstream
ss
;
ss
<<
"writing periodicity corrected mesh: "
<<
meshFile
;
print_msg
(
communicator_
,
ss
.
str
());
feMesh_
->
write_mesh
(
meshFile
);
feMesh_
->
output
(
meshFile
);
}
match
=
true
;
}
/*! \page man_mesh_write fix_modify AtC mesh write
\section syntax
fix_modify AtC mesh write <filename>
- filename = name of file to write mesh
\section examples
<TT> fix_modify AtC mesh write myMesh.mesh </TT> \n
\section description
Writes a mesh to a text file.
\section restrictions
\section related
\section default
*/
else
if
(
strcmp
(
arg
[
argIdx
],
"write"
)
==
0
)
{
argIdx
++
;
string
meshFile
=
arg
[
argIdx
];
feMesh_
->
write_mesh
(
meshFile
);
match
=
true
;
}
/*! \page man_mesh_delete_elements fix_modify AtC mesh delete_elements
\section syntax
fix_modify AtC mesh delete_elements <element_set>
- <element_set> = name of an element set
\section examples
<TT> fix_modify AtC delete_elements gap </TT>
\section description
Deletes a group of elements from the mesh.
\section restrictions
\section related
\section default
none
*/
else
if
(
strcmp
(
arg
[
argIdx
],
"delete_elements"
)
==
0
)
{
argIdx
++
;
string
esetName
=
arg
[
argIdx
];
set
<
int
>
elemSet
=
feMesh_
->
elementset
(
esetName
);
nullElements_
.
insert
(
elemSet
.
begin
(),
elemSet
.
end
());
match
=
true
;
}
else
if
(
strcmp
(
arg
[
argIdx
],
"cut"
)
==
0
)
{
argIdx
++
;
string
fsetName
=
arg
[
argIdx
++
];
set
<
PAIR
>
faceSet
=
feMesh_
->
faceset
(
fsetName
);
cutFaces_
.
insert
(
faceSet
.
begin
(),
faceSet
.
end
());
if
(
narg
>
argIdx
&&
strcmp
(
arg
[
argIdx
],
"edge"
)
==
0
)
{
argIdx
++
;
string
nsetName
=
arg
[
argIdx
];
set
<
int
>
nodeSet
=
feMesh_
->
nodeset
(
nsetName
);
cutEdge_
.
insert
(
nodeSet
.
begin
(),
nodeSet
.
end
());
}
// cut mesh
if
(
cutFaces_
.
size
()
>
0
)
cut_mesh
(
cutFaces_
,
cutEdge_
);
match
=
true
;
}
else
if
(
strcmp
(
arg
[
argIdx
],
"lammps_partition"
)
==
0
)
{
feMesh_
->
set_lammps_partition
(
true
);
match
=
true
;
}
else
if
(
strcmp
(
arg
[
argIdx
],
"data_decomposition"
)
==
0
)
{
feMesh_
->
set_data_decomposition
(
true
);
match
=
true
;
}
else
{
if
(
!
feMesh_
)
throw
ATC_Error
(
"need mesh for parsing"
);
match
=
feMesh_
->
modify
(
narg
,
arg
);
}
}
// FE_Mesh
else
{
if
(
!
feMesh_
)
throw
ATC_Error
(
"need mesh for parsing"
);
match
=
feMesh_
->
modify
(
narg
,
arg
);
}
return
match
;
}
//-----------------------------------------------------------------
void
FE_Engine
::
parse_partitions
(
int
&
argIdx
,
int
narg
,
char
**
arg
,
int
idof
,
Array
<
double
>
&
dx
)
const
{
double
x
[
3
]
=
{
0
,
0
,
0
};
int
nx
=
dx
.
size
();
// parse relative values for each element
if
(
is_numeric
(
arg
[
argIdx
]))
{
for
(
int
i
=
0
;
i
<
nx
;
++
i
)
{
if
(
is_numeric
(
arg
[
argIdx
]))
{
dx
(
i
)
=
atof
(
arg
[
argIdx
++
]);
}
else
throw
ATC_Error
(
"not enough element partitions"
);
}
}
// each segment of the piecewise funcion is length-normalized separately
else
if
(
strcmp
(
arg
[
argIdx
],
"position-number-density"
)
==
0
)
{
argIdx
++
;
double
y
[
nx
],
w
[
nx
];
int
n
[
nx
];
int
nn
=
0
;
while
(
argIdx
<
narg
)
{
if
(
!
is_numeric
(
arg
[
argIdx
]))
break
;
y
[
nn
]
=
atof
(
arg
[
argIdx
++
]);
n
[
nn
]
=
atoi
(
arg
[
argIdx
++
]);
w
[
nn
++
]
=
atof
(
arg
[
argIdx
++
]);
}
if
(
n
[
nn
-
1
]
!=
nx
)
throw
ATC_Error
(
"total element partitions do not match"
);
int
k
=
0
;
for
(
int
i
=
1
;
i
<
nn
;
++
i
)
{
int
dn
=
n
[
i
]
-
n
[
i
-
1
];
double
dy
=
y
[
i
]
-
y
[
i
-
1
];
double
w0
=
w
[
i
-
1
];
double
dw
=
w
[
i
]
-
w0
;
double
lx
=
0
;
double
l
[
dn
];
for
(
int
j
=
0
;
j
<
dn
;
++
j
)
{
double
x
=
(
j
+
0.5
)
/
dn
;
double
dl
=
w0
+
x
*
dw
;
lx
+=
dl
;
l
[
j
]
=
dl
;
}
double
scale
=
dy
/
lx
;
for
(
int
j
=
0
;
j
<
dn
;
++
j
)
{
dx
(
k
++
)
=
scale
*
l
[
j
];
}
}
}
// construct relative values from a density function
// evaluate for a domain (0,1)
else
{
XT_Function
*
f
=
XT_Function_Mgr
::
instance
()
->
function
(
&
(
arg
[
argIdx
]),
narg
-
argIdx
);
argIdx
++
;
while
(
argIdx
<
narg
)
{
if
(
!
is_numeric
(
arg
[
argIdx
++
]))
break
;
}
for
(
int
i
=
0
;
i
<
nx
;
++
i
)
{
x
[
idof
]
=
(
i
+
0.5
)
/
nx
;
dx
(
i
)
=
f
->
f
(
x
,
0.
);
}
}
}
//-----------------------------------------------------------------
void
FE_Engine
::
finish
()
{
// Nothing to do
}
//-----------------------------------------------------------------
// write geometry
//-----------------------------------------------------------------
void
FE_Engine
::
initialize_output
(
int
rank
,
string
outputPrefix
,
set
<
int
>
otypes
)
{
outputManager_
.
initialize
(
outputPrefix
,
otypes
);
if
(
!
feMesh_
)
throw
ATC_Error
(
"output needs mesh"
);
if
(
!
initialized_
)
initialize
();
if
(
!
feMesh_
->
coordinates
()
||
!
feMesh_
->
connectivity
())
throw
ATC_Error
(
"output mesh not properly initialized"
);
if
(
!
feMesh_
->
coordinates
()
->
nCols
()
||
!
feMesh_
->
connectivity
()
->
nCols
())
throw
ATC_Error
(
"output mesh is empty"
);
if
(
rank
==
0
)
outputManager_
.
write_geometry
(
feMesh_
->
coordinates
(),
feMesh_
->
connectivity
());
outputManager_
.
print_custom_names
();
}
//-----------------------------------------------------------------
// write geometry
//-----------------------------------------------------------------
void
FE_Engine
::
write_geometry
(
void
)
{
outputManager_
.
write_geometry
(
feMesh_
->
coordinates
(),
feMesh_
->
connectivity
());
}
// -------------------------------------------------------------
// write data
// -------------------------------------------------------------
void
FE_Engine
::
write_data
(
double
time
,
FIELDS
&
soln
,
OUTPUT_LIST
*
data
)
{
outputManager_
.
write_data
(
time
,
&
soln
,
data
,
(
feMesh_
->
node_map
())
->
data
());
}
// -------------------------------------------------------------
// write data
// -------------------------------------------------------------
void
FE_Engine
::
write_data
(
double
time
,
OUTPUT_LIST
*
data
)
{
outputManager_
.
write_data
(
time
,
data
,
feMesh_
->
node_map
()
->
data
());
}
// -------------------------------------------------------------
// amend mesh for deleted elements
// -------------------------------------------------------------
void
FE_Engine
::
delete_elements
(
const
set
<
int
>
&
elementList
)
{
feMesh_
->
delete_elements
(
elementList
);
}
// -------------------------------------------------------------
// amend mesh for cut at specified faces
// -------------------------------------------------------------
void
FE_Engine
::
cut_mesh
(
const
set
<
PAIR
>
&
faceSet
,
const
set
<
int
>
&
nodeSet
)
{
feMesh_
->
cut_mesh
(
faceSet
,
nodeSet
);
}
// -------------------------------------------------------------
// interpolate one value
// -------------------------------------------------------------
DENS_VEC
FE_Engine
::
interpolate_field
(
const
DENS_VEC
&
x
,
const
FIELD
&
f
)
const
{
DENS_VEC
N
;
Array
<
int
>
nodelist
;
feMesh_
->
shape_functions
(
x
,
N
,
nodelist
);
const
DENS_MAT
&
vI
=
f
.
quantity
();
int
dof
=
vI
.
nCols
();
DENS_MAT
vIe
(
nNodesPerElement_
,
dof
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
i
++
)
for
(
int
j
=
0
;
j
<
dof
;
j
++
)
vIe
(
i
,
j
)
=
vI
(
nodelist
(
i
),
j
);
DENS_VEC
vP
;
vP
=
N
*
vIe
;
return
vP
;
}
// -------------------------------------------------------------
// interpolate fields and gradients
// Currently, this function will break if called with an unowned ielem.
// Currently, this function is only called with owned ielems.
// -------------------------------------------------------------
void
FE_Engine
::
interpolate_fields
(
const
int
ielem
,
const
FIELDS
&
fields
,
AliasArray
<
int
>
&
conn
,
DENS_MAT
&
N
,
DENS_MAT_VEC
&
dN
,
DIAG_MAT
&
_weights_
,
FIELD_MATS
&
fieldsAtIPs
,
GRAD_FIELD_MATS
&
gradFieldsAtIPs
)
const
{
// evaluate shape functions
feMesh_
->
shape_function
(
ielem
,
N
,
dN
,
_weights_
);
// get connectivity
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
// compute fields and gradients of fields ips of this element
FIELD_MATS
localElementFields
;
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
// field values at all nodes
_fieldName_
=
_fieldItr_
->
first
;
const
DENS_MAT
&
vI
=
(
_fieldItr_
->
second
).
quantity
();
int
dof
=
vI
.
nCols
();
// field values at integration points -> to be computed
DENS_MAT
&
vP
=
fieldsAtIPs
[
_fieldName_
];
// gradients of field at integration points -> to be computed
DENS_MAT_VEC
&
dvP
=
gradFieldsAtIPs
[
_fieldName_
];
if
(
_fieldName_
==
ELECTRON_WAVEFUNCTION_ENERGIES
)
{
vP
=
vI
;
continue
;
}
// field values at element nodes
DENS_MAT
&
vIe
=
localElementFields
[
_fieldName_
];
// gather local field
vIe
.
reset
(
nNodesPerElement_
,
dof
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
i
++
)
for
(
int
j
=
0
;
j
<
dof
;
j
++
)
vIe
(
i
,
j
)
=
vI
(
conn
(
i
),
j
);
// interpolate field at integration points
vP
=
N
*
vIe
;
// gradients
dvP
.
assign
(
nSD_
,
DENS_MAT
(
nIPsPerElement_
,
dof
));
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
dvP
[
j
]
=
dN
[
j
]
*
vIe
;
}
}
// -------------------------------------------------------------
// interpolate fields
// Currently, this function will break if called with an unowned ielem.
// Currently, this function is only called with owned ielems.
// -------------------------------------------------------------
void
FE_Engine
::
interpolate_fields
(
const
int
ielem
,
const
FIELDS
&
fields
,
AliasArray
<
int
>
&
conn
,
DENS_MAT
&
N
,
DIAG_MAT
&
_weights_
,
FIELD_MATS
&
fieldsAtIPs
)
const
{
// evaluate shape functions
feMesh_
->
shape_function
(
ielem
,
N
,
_weights_
);
// get connectivity
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
// compute fields and gradients of fields ips of this element
FIELD_MATS
localElementFields
;
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
// field values at all nodes
_fieldName_
=
_fieldItr_
->
first
;
const
DENS_MAT
&
vI
=
(
_fieldItr_
->
second
).
quantity
();
int
dof
=
vI
.
nCols
();
// field values at integration points -> to be computed
DENS_MAT
&
vP
=
fieldsAtIPs
[
_fieldName_
];
// field values at element nodes
DENS_MAT
&
vIe
=
localElementFields
[
_fieldName_
];
if
(
_fieldName_
==
ELECTRON_WAVEFUNCTION_ENERGIES
)
{
vP
=
vI
;
continue
;
}
// gather local field
vIe
.
reset
(
nNodesPerElement_
,
dof
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
i
++
)
for
(
int
j
=
0
;
j
<
dof
;
j
++
)
vIe
(
i
,
j
)
=
vI
(
conn
(
i
),
j
);
// interpolate field at integration points
vP
=
N
*
vIe
;
}
}
// -------------------------------------------------------------
// compute dimensionless stiffness matrix using native quadrature
// -------------------------------------------------------------
void
FE_Engine
::
stiffness_matrix
(
SPAR_MAT
&
matrix
)
const
{
// assemble consistent mass
matrix
.
reset
(
nNodesUnique_
,
nNodesUnique_
);
// zero since partial fill
DENS_MAT
elementMassMatrix
(
nNodesPerElement_
,
nNodesPerElement_
);
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// evaluate shape functions
feMesh_
->
shape_function
(
ielem
,
_N_
,
_dN_
,
_weights_
);
// _N_ unused
// perform quadrature
elementMassMatrix
=
_dN_
[
0
].
transMat
(
_weights_
*
_dN_
[
0
]);
for
(
int
i
=
1
;
i
<
nSD_
;
++
i
)
{
elementMassMatrix
+=
_dN_
[
i
].
transMat
(
_weights_
*
_dN_
[
i
]);
}
// get connectivity
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
_conn_
(
j
);
matrix
.
add
(
inode
,
jnode
,
elementMassMatrix
(
i
,
j
));
}
}
}
#ifdef ISOLATE_FE
sparse_allsum
(
communicator_
,
matrix
);
#else
LammpsInterface
::
instance
()
->
sparse_allsum
(
matrix
);
#endif
matrix
.
compress
();
}
// -------------------------------------------------------------
// compute tangent using native quadrature for one (field,field) pair
// -------------------------------------------------------------
void
FE_Engine
::
compute_tangent_matrix
(
const
Array2D
<
bool
>
&
rhsMask
,
const
pair
<
FieldName
,
FieldName
>
row_col
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
int
>
&
elementMaterials
,
SPAR_MAT
&
tangent
,
const
DenseMatrix
<
bool
>
*
elementMask
)
const
{
tangent
.
reset
(
nNodesUnique_
,
nNodesUnique_
);
FieldName
rowField
=
row_col
.
first
;
FieldName
colField
=
row_col
.
second
;
bool
BB
=
rhsMask
(
rowField
,
FLUX
);
bool
NN
=
rhsMask
(
rowField
,
SOURCE
);
DENS_MAT
elementMatrix
(
nNodesPerElement_
,
nNodesPerElement_
);
DENS_MAT
coefsAtIPs
;
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// if element is masked, skip it
if
(
elementMask
&&
!
(
*
elementMask
)(
ielem
,
0
))
continue
;
// material id
int
imat
=
elementMaterials
(
ielem
);
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
// interpolate fields and gradients (nonlinear only)
interpolate_fields
(
ielem
,
fields
,
_conn_
,
_N_
,
_dN_
,
_weights_
,
_fieldsAtIPs_
,
_gradFieldsAtIPs_
);
// evaluate Physics model
if
(
!
(
physicsModel
->
null
(
rowField
,
imat
))
)
{
if
(
BB
&&
physicsModel
->
weak_equation
(
rowField
)
->
has_BB_tangent_coefficients
()
)
{
physicsModel
->
weak_equation
(
rowField
)
->
BB_tangent_coefficients
(
colField
,
_fieldsAtIPs_
,
mat
,
coefsAtIPs
);
DIAG_MAT
D
(
column
(
coefsAtIPs
,
0
));
D
=
_weights_
*
D
;
elementMatrix
=
_dN_
[
0
].
transMat
(
D
*
_dN_
[
0
]);
for
(
int
i
=
1
;
i
<
nSD_
;
i
++
)
{
elementMatrix
+=
_dN_
[
i
].
transMat
(
D
*
_dN_
[
i
]);
}
}
else
{
elementMatrix
.
reset
(
nNodesPerElement_
,
nNodesPerElement_
);
}
if
(
NN
&&
physicsModel
->
weak_equation
(
rowField
)
->
has_NN_tangent_coefficients
()
)
{
physicsModel
->
weak_equation
(
rowField
)
->
NN_tangent_coefficients
(
colField
,
_fieldsAtIPs_
,
mat
,
coefsAtIPs
);
DIAG_MAT
D
(
column
(
coefsAtIPs
,
0
));
D
=
_weights_
*
D
;
elementMatrix
+=
_N_
.
transMat
(
D
*
_N_
);
}
// assemble
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
_conn_
(
j
);
tangent
.
add
(
inode
,
jnode
,
elementMatrix
(
i
,
j
));
}
}
}
}
#ifdef ISOLATE_FE
sparse_allsum
(
communicator_
,
tangent
);
#else
LammpsInterface
::
instance
()
->
sparse_allsum
(
tangent
);
#endif
tangent
.
compress
();
}
// -------------------------------------------------------------
// compute tangent using given quadrature for one (field,field) pair
// -------------------------------------------------------------
void
FE_Engine
::
compute_tangent_matrix
(
const
Array2D
<
bool
>
&
rhsMask
,
const
pair
<
FieldName
,
FieldName
>
row_col
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
set
<
int
>
>
&
pointMaterialGroups
,
const
DIAG_MAT
&
weights
,
const
SPAR_MAT
&
N
,
const
SPAR_MAT_VEC
&
dN
,
SPAR_MAT
&
tangent
,
const
DenseMatrix
<
bool
>
*
elementMask
)
const
{
int
nn
=
nNodesUnique_
;
FieldName
rowField
=
row_col
.
first
;
FieldName
colField
=
row_col
.
second
;
bool
BB
=
rhsMask
(
rowField
,
FLUX
);
bool
NN
=
rhsMask
(
rowField
,
SOURCE
);
DENS_MAT
K
(
nn
,
nn
);
DENS_MAT
coefsAtIPs
;
int
nips
=
weights
.
nCols
();
if
(
nips
>
0
)
{
// compute fields and gradients of fields at given ips
GRAD_FIELD_MATS
gradFieldsAtIPs
;
FIELD_MATS
fieldsAtIPs
;
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
_fieldName_
=
_fieldItr_
->
first
;
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
dof
=
field
.
nCols
();
gradFieldsAtIPs
[
_fieldName_
].
assign
(
nSD_
,
DENS_MAT
(
nips
,
dof
));
fieldsAtIPs
[
_fieldName_
]
=
N
*
field
;
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
gradFieldsAtIPs
[
_fieldName_
][
j
]
=
(
*
dN
[
j
])
*
field
;
}
}
// treat single material point sets specially
int
nMatls
=
pointMaterialGroups
.
size
();
int
atomMatls
=
0
;
for
(
int
imat
=
0
;
imat
<
nMatls
;
imat
++
)
{
const
set
<
int
>
&
indices
=
pointMaterialGroups
(
imat
);
if
(
indices
.
size
()
>
0
)
atomMatls
++
;
}
bool
singleMaterial
=
(
atomMatls
==
1
);
if
(
!
singleMaterial
)
throw
ATC_Error
(
"FE_Engine::compute_tangent_matrix-given quadrature can not handle multiple atom material currently"
);
if
(
singleMaterial
)
{
int
imat
=
0
;
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
// evaluate Physics model
if
(
!
(
physicsModel
->
null
(
rowField
,
imat
))
)
{
if
(
BB
&&
physicsModel
->
weak_equation
(
rowField
)
->
has_BB_tangent_coefficients
()
)
{
physicsModel
->
weak_equation
(
rowField
)
->
BB_tangent_coefficients
(
colField
,
fieldsAtIPs
,
mat
,
coefsAtIPs
);
DIAG_MAT
D
(
column
(
coefsAtIPs
,
0
));
D
=
weights
*
D
;
K
=
(
*
dN
[
0
]).
transMat
(
D
*
(
*
dN
[
0
]));
for
(
int
i
=
1
;
i
<
nSD_
;
i
++
)
{
K
+=
(
*
dN
[
i
]).
transMat
(
D
*
(
*
dN
[
i
]));
}
}
if
(
NN
&&
physicsModel
->
weak_equation
(
rowField
)
->
has_NN_tangent_coefficients
()
)
{
physicsModel
->
weak_equation
(
rowField
)
->
NN_tangent_coefficients
(
colField
,
fieldsAtIPs
,
mat
,
coefsAtIPs
);
DIAG_MAT
D
(
column
(
coefsAtIPs
,
0
));
D
=
weights
*
D
;
K
+=
N
.
transMat
(
D
*
N
);
}
}
}
}
// share information between processors
int
count
=
nn
*
nn
;
DENS_MAT
K_sum
(
nn
,
nn
);
allsum
(
communicator_
,
K
.
ptr
(),
K_sum
.
ptr
(),
count
);
// create sparse from dense
tangent
.
reset
(
K_sum
,
tol_sparse
);
tangent
.
compress
();
}
// -------------------------------------------------------------
// compute a consistent mass matrix for native quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_mass_matrix
(
const
Array
<
FieldName
>&
field_mask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
int
>
&
elementMaterials
,
CON_MASS_MATS
&
massMatrices
,
const
DenseMatrix
<
bool
>
*
elementMask
)
const
{
int
nfields
=
field_mask
.
size
();
vector
<
FieldName
>
usedFields
;
DENS_MAT
elementMassMatrix
(
nNodesPerElement_
,
nNodesPerElement_
);
// (JAT, 04/21/11) FIX THIS
DENS_MAT
capacity
;
// zero, use incoming matrix as template if possible
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
SPAR_MAT
&
M
=
massMatrices
[
_fieldName_
].
set_quantity
();
// compresses 2May11
if
(
M
.
has_template
())
{
M
=
0
;
}
else
{
M
.
reset
(
nNodesUnique_
,
nNodesUnique_
);
}
M
.
reset
(
nNodesUnique_
,
nNodesUnique_
);
}
// element wise assembly
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// if element is masked, skip
if
(
elementMask
&&
!
(
*
elementMask
)(
ielem
,
0
))
continue
;
// material id
int
imat
=
elementMaterials
(
ielem
);
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
// interpolate fields
interpolate_fields
(
ielem
,
fields
,
_conn_
,
_N_
,
_weights_
,
_fieldsAtIPs_
);
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
SPAR_MAT
&
M
=
massMatrices
[
_fieldName_
].
set_quantity
();
// skip null weakEqns by material
if
(
!
(
physicsModel
->
null
(
_fieldName_
,
imat
))
)
{
physicsModel
->
weak_equation
(
_fieldName_
)
->
M_integrand
(
_fieldsAtIPs_
,
mat
,
capacity
);
DIAG_MAT
rho
(
column
(
capacity
,
0
));
elementMassMatrix
=
_N_
.
transMat
(
_weights_
*
rho
*
_N_
);
// assemble
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
_conn_
(
j
);
M
.
add
(
inode
,
jnode
,
elementMassMatrix
(
i
,
j
));
}
}
}
}
}
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
SPAR_MAT
&
M
=
massMatrices
[
_fieldName_
].
set_quantity
();
#ifdef ISOLATE_FE
sparse_allsum
(
communicator_
,
M
);
#else
LammpsInterface
::
instance
()
->
sparse_allsum
(
M
);
#endif
}
// fix zero diagonal entries due to null material elements
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
SPAR_MAT
&
M
=
massMatrices
[
_fieldName_
].
set_quantity
();
for
(
int
inode
=
0
;
inode
<
nNodesUnique_
;
++
inode
)
{
if
(
!
M
.
has_entry
(
inode
,
inode
))
{
M
.
set
(
inode
,
inode
,
1.0
);
}
}
M
.
compress
();
}
}
// -------------------------------------------------------------
// compute dimensionless consistent mass using native quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_mass_matrix
(
SPAR_MAT
&
massMatrix
)
const
{
// assemble nnodes X nnodes matrix
massMatrix
.
reset
(
nNodesUnique_
,
nNodesUnique_
);
// zero since partial fill
DENS_MAT
elementMassMatrix
(
nNodesPerElement_
,
nNodesPerElement_
);
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// evaluate shape functions
feMesh_
->
shape_function
(
ielem
,
_N_
,
_weights_
);
// perform quadrature
elementMassMatrix
=
_N_
.
transMat
(
_weights_
*
_N_
);
// get connectivity
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
_conn_
(
j
);
massMatrix
.
add
(
inode
,
jnode
,
elementMassMatrix
(
i
,
j
));
}
}
}
// Assemble partial results
#ifdef ISOLATE_FE
sparse_allsum
(
communicator_
,
massMatrix
);
#else
LammpsInterface
::
instance
()
->
sparse_allsum
(
massMatrix
);
#endif
}
// -------------------------------------------------------------
// compute dimensionless consistent mass using given quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_mass_matrix
(
const
DIAG_MAT
&
weights
,
const
SPAR_MAT
&
N
,
SPAR_MAT
&
massMatrix
)
const
{
int
nn
=
N
.
nCols
();
int
nips
=
N
.
nRows
();
DENS_MAT
tmp_mass_matrix_local
(
nn
,
nn
),
tmp_mass_matrix
(
nn
,
nn
);
if
(
nips
>
0
)
{
tmp_mass_matrix_local
=
N
.
transMat
(
weights
*
N
);
}
// share information between processors
int
count
=
nn
*
nn
;
allsum
(
communicator_
,
tmp_mass_matrix_local
.
ptr
(),
tmp_mass_matrix
.
ptr
(),
count
);
// create sparse from dense
massMatrix
.
reset
(
tmp_mass_matrix
,
tol_sparse
);
}
// -------------------------------------------------------------
// compute dimensionless lumped mass using native quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_lumped_mass_matrix
(
DIAG_MAT
&
M
)
const
{
M
.
reset
(
nNodesUnique_
,
nNodesUnique_
);
// zero since partial fill
// assemble lumped diagonal mass
DENS_VEC
Nvec
(
nNodesPerElement_
);
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// evaluate shape functions
feMesh_
->
shape_function
(
ielem
,
_N_
,
_weights_
);
CLON_VEC
w
(
_weights_
);
// get connectivity
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
Nvec
=
w
*
_N_
;
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
M
(
inode
,
inode
)
+=
Nvec
(
i
);
}
}
// Assemble partial results
allsum
(
communicator_
,
MPI_IN_PLACE
,
M
.
ptr
(),
M
.
size
());
}
// -------------------------------------------------------------
// compute physical lumped mass using native quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_lumped_mass_matrix
(
const
Array
<
FieldName
>&
field_mask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
int
>
&
elementMaterials
,
MASS_MATS
&
massMatrices
,
// mass matrix
const
DenseMatrix
<
bool
>
*
elementMask
)
const
{
int
nfields
=
field_mask
.
size
();
// zero initialize for assembly
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
DIAG_MAT
&
M
=
massMatrices
[
field_mask
(
j
)].
set_quantity
();
M
.
reset
(
nNodesUnique_
,
nNodesUnique_
);
}
// assemble diagonal matrix
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// if element is masked, skip it
if
(
elementMask
&&
!
(
*
elementMask
)(
ielem
,
0
))
continue
;
// material id
int
imat
=
elementMaterials
(
ielem
);
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
interpolate_fields
(
ielem
,
fields
,
_conn_
,
_N_
,
_weights_
,
_fieldsAtIPs_
);
// compute densities, integrate & assemble
DENS_MAT
capacity
;
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
if
(
!
physicsModel
->
null
(
_fieldName_
,
imat
))
{
physicsModel
->
weak_equation
(
_fieldName_
)
->
M_integrand
(
_fieldsAtIPs_
,
mat
,
capacity
);
_Nmat_
=
_N_
.
transMat
(
_weights_
*
capacity
);
DIAG_MAT
&
M
=
massMatrices
[
_fieldName_
].
set_quantity
();
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
M
(
inode
,
inode
)
+=
_Nmat_
(
i
,
0
);
}
}
}
}
// Assemble partial results
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
DIAG_MAT
&
M
=
massMatrices
[
_fieldName_
].
set_quantity
();
allsum
(
communicator_
,
MPI_IN_PLACE
,
M
.
ptr
(),
M
.
size
());
}
// fix zero diagonal entries due to null material elements
for
(
int
inode
=
0
;
inode
<
nNodesUnique_
;
++
inode
)
{
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
DIAG_MAT
&
M
=
massMatrices
[
_fieldName_
].
set_quantity
();
if
(
M
(
inode
,
inode
)
==
0.0
)
{
M
(
inode
,
inode
)
=
1.0
;
}
}
}
}
// -------------------------------------------------------------
// compute physical lumped mass using given quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_lumped_mass_matrix
(
const
Array
<
FieldName
>
&
field_mask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
set
<
int
>
>
&
pointMaterialGroups
,
const
DIAG_MAT
&
weights
,
const
SPAR_MAT
&
N
,
MASS_MATS
&
M
)
const
// mass matrices
{
int
nips
=
weights
.
nCols
();
int
nfields
=
field_mask
.
size
();
// initialize
map
<
FieldName
,
DIAG_MAT
>
M_local
;
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
M_local
[
_fieldName_
].
reset
(
nNodesUnique_
,
nNodesUnique_
);
M
[
_fieldName_
].
reset
(
nNodesUnique_
,
nNodesUnique_
);
}
if
(
nips
>
0
)
{
// compute fields at given ips
// compute all fields since we don't the capacities dependencies
FIELD_MATS
fieldsAtIPs
;
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
_fieldName_
=
_fieldItr_
->
first
;
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
dof
=
field
.
nCols
();
fieldsAtIPs
[
_fieldName_
].
reset
(
nips
,
dof
);
fieldsAtIPs
[
_fieldName_
]
=
N
*
field
;
}
// treat single material point sets specially
int
nMatls
=
pointMaterialGroups
.
size
();
int
atomMatls
=
0
;
int
iatomMat
=
0
;
for
(
int
imat
=
0
;
imat
<
nMatls
;
imat
++
)
{
const
set
<
int
>
&
indices
=
pointMaterialGroups
(
imat
);
if
(
indices
.
size
()
>
0
)
{
iatomMat
=
imat
;
atomMatls
++
;
}
}
if
(
atomMatls
==
0
)
{
throw
ATC_Error
(
"no materials in atom region - atoms may have migrated to FE-only region"
);
}
bool
singleMaterial
=
(
atomMatls
==
1
);
if
(
!
singleMaterial
)
{
stringstream
ss
;
ss
<<
" WARNING: multiple materials in atomic region"
;
print_msg
(
communicator_
,
ss
.
str
());
}
// setup data structures
FIELD_MATS
capacities
;
// evaluate physics model & compute capacity|density for requested fields
if
(
singleMaterial
)
{
const
Material
*
mat
=
physicsModel
->
material
(
iatomMat
);
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
// mass matrix needs to be invertible so null matls have cap=1
if
(
physicsModel
->
null
(
_fieldName_
,
iatomMat
))
{
throw
ATC_Error
(
"null material not supported for atomic region (single material)"
);
const
FIELD
&
f
=
(
fields
.
find
(
_fieldName_
))
->
second
;
capacities
[
_fieldName_
].
reset
(
f
.
nRows
(),
f
.
nCols
());
capacities
[
_fieldName_
]
=
1.
;
}
else
{
physicsModel
->
weak_equation
(
_fieldName_
)
->
M_integrand
(
fieldsAtIPs
,
mat
,
capacities
[
_fieldName_
]);
}
}
}
else
{
FIELD_MATS
groupCapacities
,
fieldsGroup
;
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
capacities
[
_fieldName_
].
reset
(
nips
,
1
);
}
for
(
int
imat
=
0
;
imat
<
pointMaterialGroups
.
size
();
imat
++
)
{
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
const
set
<
int
>
&
indices
=
pointMaterialGroups
(
imat
);
int
npts
=
indices
.
size
();
if
(
npts
>
0
)
{
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
groupCapacities
[
_fieldName_
].
reset
(
npts
,
1
);
const
FIELD
&
f
=
(
fields
.
find
(
_fieldName_
))
->
second
;
int
ndof
=
f
.
nCols
();
fieldsGroup
[
_fieldName_
].
reset
(
npts
,
ndof
);
int
i
=
0
;
for
(
set
<
int
>::
const_iterator
iter
=
indices
.
begin
();
iter
!=
indices
.
end
();
iter
++
,
i
++
)
{
for
(
int
dof
=
0
;
dof
<
ndof
;
++
dof
)
{
fieldsGroup
[
_fieldName_
](
i
,
dof
)
=
fieldsAtIPs
[
_fieldName_
](
*
iter
,
dof
);
}
}
}
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
if
(
physicsModel
->
null
(
_fieldName_
,
imat
))
{
throw
ATC_Error
(
"null material not supported for atomic region (multiple materials)"
);
const
FIELD
&
f
=
(
fields
.
find
(
_fieldName_
))
->
second
;
groupCapacities
[
_fieldName_
].
reset
(
f
.
nRows
(),
f
.
nCols
());
groupCapacities
[
_fieldName_
]
=
1.
;
}
else
{
physicsModel
->
weak_equation
(
_fieldName_
)
->
M_integrand
(
fieldsGroup
,
mat
,
groupCapacities
[
_fieldName_
]);
}
int
i
=
0
;
for
(
set
<
int
>::
const_iterator
iter
=
indices
.
begin
();
iter
!=
indices
.
end
();
iter
++
,
i
++
)
{
capacities
[
_fieldName_
](
*
iter
,
0
)
=
groupCapacities
[
_fieldName_
](
i
,
0
);
}
}
}
}
}
// integrate & assemble
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
M_local
[
_fieldName_
].
reset
(
// assume all columns same
column
(
N
.
transMat
(
weights
*
capacities
[
_fieldName_
]),
0
)
);
}
}
// Share information between processors
for
(
int
j
=
0
;
j
<
nfields
;
++
j
)
{
_fieldName_
=
field_mask
(
j
);
DIAG_MAT
&
myMassMat
(
M
[
_fieldName_
].
set_quantity
());
int
count
=
M_local
[
_fieldName_
].
size
();
allsum
(
communicator_
,
M_local
[
_fieldName_
].
ptr
(),
myMassMat
.
ptr
(),
count
);
}
}
//-----------------------------------------------------------------
// compute assembled average gradient evaluated at the nodes
//-----------------------------------------------------------------
void
FE_Engine
::
compute_gradient_matrix
(
SPAR_MAT_VEC
&
grad_matrix
)
const
{
// zero
DENS_MAT_VEC
tmp_grad_matrix
(
nSD_
);
for
(
int
i
=
0
;
i
<
nSD_
;
i
++
)
{
tmp_grad_matrix
[
i
].
reset
(
nNodesUnique_
,
nNodesUnique_
);
}
// element wise assembly
Array
<
int
>
count
(
nNodesUnique_
);
count
=
0
;
set_quadrature
(
NODAL
);
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// evaluate shape functions
feMesh_
->
shape_function
(
ielem
,
_N_
,
_dN_
,
_weights_
);
// get connectivity
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
for
(
int
j
=
0
;
j
<
nIPsPerElement_
;
++
j
)
{
int
jnode
=
_conn_
(
j
);
count
(
jnode
)
+=
1
;
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
tmp_grad_matrix
[
k
](
jnode
,
inode
)
+=
_dN_
[
k
](
j
,
i
);
}
}
}
}
// Assemble partial results
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
allsum
(
communicator_
,
MPI_IN_PLACE
,
tmp_grad_matrix
[
k
].
ptr
(),
tmp_grad_matrix
[
k
].
size
());
}
int_allsum
(
communicator_
,
MPI_IN_PLACE
,
count
.
ptr
(),
count
.
size
());
set_quadrature
(
quadrature_
);
//reset to default
for
(
int
inode
=
0
;
inode
<
nNodesUnique_
;
++
inode
)
{
for
(
int
jnode
=
0
;
jnode
<
nNodesUnique_
;
++
jnode
)
{
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
tmp_grad_matrix
[
k
](
jnode
,
inode
)
/=
count
(
jnode
);
}
}
}
// compact dense matrices
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
grad_matrix
[
k
]
->
reset
(
tmp_grad_matrix
[
k
],
tol_sparse
);
}
}
// -------------------------------------------------------------
// compute energy per node using native quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_energy
(
const
Array
<
FieldName
>
&
mask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
int
>
&
elementMaterials
,
FIELD_MATS
&
energies
,
const
DenseMatrix
<
bool
>
*
elementMask
,
const
IntegrationDomainType
domain
)
const
{
// Zero out all fields
for
(
int
n
=
0
;
n
<
mask
.
size
();
n
++
)
{
_fieldName_
=
mask
(
n
);
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
energies
[
_fieldName_
].
reset
(
field
.
nRows
(),
1
);
}
DENS_MAT
elementEnergy
(
nNodesPerElement_
,
1
);
// workspace
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// if element is masked, skip it
if
(
domain
!=
FULL_DOMAIN
&&
elementMask
&&
!
(
*
elementMask
)(
ielem
,
0
))
continue
;
// material id
int
imat
=
elementMaterials
(
ielem
);
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
interpolate_fields
(
ielem
,
fields
,
_conn_
,
_N_
,
_dN_
,
_weights_
,
_fieldsAtIPs_
,
_gradFieldsAtIPs_
);
// assemble
for
(
int
n
=
0
;
n
<
mask
.
size
();
n
++
)
{
_fieldName_
=
mask
(
n
);
if
(
physicsModel
->
null
(
_fieldName_
,
imat
))
continue
;
if
(
!
(
physicsModel
->
weak_equation
(
_fieldName_
)
->
has_E_integrand
()))
continue
;
physicsModel
->
weak_equation
(
_fieldName_
)
->
E_integrand
(
_fieldsAtIPs_
,
_gradFieldsAtIPs_
,
mat
,
elementEnergy
);
_fieldItr_
=
fields
.
find
(
_fieldName_
);
_Nmat_
=
_N_
.
transMat
(
_weights_
*
elementEnergy
);
DENS_MAT
&
energy
=
energies
[
_fieldName_
];
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
energy
(
inode
,
0
)
+=
_Nmat_
(
i
,
0
);
}
}
}
for
(
int
n
=
0
;
n
<
mask
.
size
();
n
++
)
{
_fieldName_
=
mask
(
n
);
DENS_MAT
&
myEnergy
(
energies
[
_fieldName_
]);
allsum
(
communicator_
,
MPI_IN_PLACE
,
myEnergy
.
ptr
(),
myEnergy
.
size
());
}
}
// -------------------------------------------------------------
// compute rhs using native quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_rhs_vector
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
int
>
&
elementMaterials
,
FIELDS
&
rhs
,
bool
freeOnly
,
const
DenseMatrix
<
bool
>
*
elementMask
)
const
{
vector
<
FieldName
>
usedFields
;
// size and zero output
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
rhsMask
(
_fieldName_
,
FLUX
)
||
rhsMask
(
_fieldName_
,
SOURCE
))
{
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
rhs
[
_fieldName_
].
reset
(
field
.
nRows
(),
field
.
nCols
());
// Save field names for easy lookup later.
usedFields
.
push_back
(
_fieldName_
);
}
}
// Iterate over elements partitioned onto current processor.
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// Skip masked elements
if
(
elementMask
&&
!
(
*
elementMask
)(
ielem
,
0
))
continue
;
int
imat
=
elementMaterials
(
ielem
);
// material id
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
// interpolate field values to integration points
interpolate_fields
(
ielem
,
fields
,
_conn_
,
_N_
,
_dN_
,
_weights_
,
_fieldsAtIPs_
,
_gradFieldsAtIPs_
);
// rescale by _weights_, a diagonal matrix
_Nw_
=
_weights_
*
_N_
;
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
_dNw_
[
j
]
=
_weights_
*
_dN_
[
j
];
// evaluate physics model and assemble
// _Nfluxes is a set of fields
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
!
rhsMask
(
_fieldName_
,
SOURCE
))
continue
;
if
(
physicsModel
->
null
(
_fieldName_
,
imat
))
continue
;
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
DENS_MAT
&
myRhs
(
rhs
[
_fieldName_
].
set_quantity
());
int
dof
=
field
.
nCols
();
bool
has
=
physicsModel
->
weak_equation
(
_fieldName_
)
->
N_integrand
(
_fieldsAtIPs_
,
_gradFieldsAtIPs_
,
mat
,
_Nfluxes_
);
if
(
!
has
)
continue
;
_Nmat_
=
_Nw_
.
transMat
(
_Nfluxes_
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
k
=
0
;
k
<
dof
;
++
k
)
{
myRhs
(
inode
,
k
)
+=
_Nmat_
(
i
,
k
);
}
}
}
// _Bfluxes_ is a set of field gradients
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
!
rhsMask
(
_fieldName_
,
FLUX
))
continue
;
if
(
physicsModel
->
null
(
_fieldName_
,
imat
))
continue
;
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
DENS_MAT
&
myRhs
(
rhs
[
_fieldName_
].
set_quantity
());
int
dof
=
field
.
nCols
();
physicsModel
->
weak_equation
(
_fieldName_
)
->
B_integrand
(
_fieldsAtIPs_
,
_gradFieldsAtIPs_
,
mat
,
_Bfluxes_
);
for
(
int
j
=
0
;
j
<
nSD_
;
j
++
)
{
_Nmat_
=
_dNw_
[
j
].
transMat
(
_Bfluxes_
[
j
]);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
k
=
0
;
k
<
dof
;
++
k
)
{
myRhs
(
inode
,
k
)
+=
_Nmat_
(
i
,
k
);
}
}
}
}
}
vector
<
FieldName
>::
iterator
_fieldIter_
;
for
(
_fieldIter_
=
usedFields
.
begin
();
_fieldIter_
!=
usedFields
.
end
();
++
_fieldIter_
)
{
// myRhs is where we put the global result for this field.
_fieldName_
=
*
_fieldIter_
;
DENS_MAT
&
myRhs
(
rhs
[
_fieldName_
].
set_quantity
());
// Sum matrices in-place across all processors into myRhs.
allsum
(
communicator_
,
MPI_IN_PLACE
,
myRhs
.
ptr
(),
myRhs
.
size
());
}
}
// -------------------------------------------------------------
// compute rhs using given quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_rhs_vector
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
set
<
int
>
>&
pointMaterialGroups
,
const
DIAG_MAT
&
weights
,
const
SPAR_MAT
&
N
,
const
SPAR_MAT_VEC
&
dN
,
FIELDS
&
rhs
)
const
{
FIELD_MATS
rhs_local
;
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
rhsMask
(
_fieldName_
,
FLUX
)
||
rhsMask
(
_fieldName_
,
SOURCE
))
{
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
nrows
=
field
.
nRows
();
int
ncols
=
field
.
nCols
();
rhs
[
_fieldName_
].
reset
(
nrows
,
ncols
);
rhs_local
[
_fieldName_
].
reset
(
nrows
,
ncols
);
}
}
int
nips
=
weights
.
nCols
();
if
(
nips
>
0
)
{
// compute fields and gradients of fields at given ips
GRAD_FIELD_MATS
gradFieldsAtIPs
;
FIELD_MATS
fieldsAtIPs
;
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
_fieldName_
=
_fieldItr_
->
first
;
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
dof
=
field
.
nCols
();
gradFieldsAtIPs
[
_fieldName_
].
assign
(
nSD_
,
DENS_MAT
(
nips
,
dof
));
fieldsAtIPs
[
_fieldName_
]
=
N
*
field
;
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
gradFieldsAtIPs
[
_fieldName_
][
j
]
=
(
*
dN
[
j
])
*
field
;
}
}
// setup data structures
FIELD_MATS
Nfluxes
;
GRAD_FIELD_MATS
Bfluxes
;
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
int
index
=
(
int
)
_fieldName_
;
if
(
rhsMask
(
index
,
FLUX
)
)
{
Bfluxes
[
_fieldName_
].
assign
(
nSD_
,
DENS_MAT
());
}
}
// treat single material point sets specially
int
nMatls
=
pointMaterialGroups
.
size
();
int
atomMatls
=
0
;
for
(
int
imat
=
0
;
imat
<
nMatls
;
imat
++
)
{
const
set
<
int
>
&
indices
=
pointMaterialGroups
(
imat
);
if
(
indices
.
size
()
>
0
)
atomMatls
++
;
}
bool
singleMaterial
=
(
atomMatls
==
1
);
// evaluate physics model
if
(
singleMaterial
)
{
int
imat
=
0
;
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
int
index
=
(
int
)
_fieldName_
;
if
(
!
physicsModel
->
null
(
_fieldName_
,
imat
))
{
if
(
rhsMask
(
index
,
SOURCE
)
)
{
physicsModel
->
weak_equation
(
_fieldName_
)
->
N_integrand
(
fieldsAtIPs
,
gradFieldsAtIPs
,
mat
,
Nfluxes
[
_fieldName_
]);
}
if
(
rhsMask
(
index
,
FLUX
)
)
{
physicsModel
->
weak_equation
(
_fieldName_
)
->
B_integrand
(
fieldsAtIPs
,
gradFieldsAtIPs
,
mat
,
Bfluxes
[
_fieldName_
]);
}
}
}
}
else
{
// 1) permanent workspace with per-row mapped clones per matl
// from caller/atc
// 2) per point calls to N/B_integrand
// 3) collect internalToAtom into materials and use mem-cont clones
// what about dof for matrices and data storage: clone _matrix_
// for each material group:
// set up storage
DENS_MAT
group_Nfluxes
;
DENS_MAT_VEC
group_Bfluxes
;
group_Bfluxes
.
reserve
(
nSD_
);
FIELD_MATS
fieldsGroup
;
GRAD_FIELD_MATS
gradFieldsGroup
;
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
_fieldName_
=
_fieldItr_
->
first
;
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
ndof
=
field
.
nCols
();
gradFieldsGroup
[
_fieldName_
].
assign
(
nSD_
,
DENS_MAT
(
nips
,
ndof
));
Nfluxes
[
_fieldName_
].
reset
(
nips
,
ndof
);
Bfluxes
[
_fieldName_
].
assign
(
nSD_
,
DENS_MAT
(
nips
,
ndof
));
//}
}
// copy fields
for
(
int
imat
=
0
;
imat
<
pointMaterialGroups
.
size
();
imat
++
)
{
const
set
<
int
>
&
indices
=
pointMaterialGroups
(
imat
);
const
Material
*
mat
=
physicsModel
->
material
(
0
);
int
npts
=
indices
.
size
();
int
i
=
0
;
for
(
set
<
int
>::
const_iterator
iter
=
indices
.
begin
();
iter
!=
indices
.
end
();
iter
++
,
i
++
)
{
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
_fieldName_
=
_fieldItr_
->
first
;
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
ndof
=
field
.
nCols
();
fieldsGroup
[
_fieldName_
].
reset
(
npts
,
ndof
);
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
(
gradFieldsGroup
[
_fieldName_
][
j
]).
reset
(
npts
,
ndof
);
}
for
(
int
dof
=
0
;
dof
<
ndof
;
++
dof
)
{
fieldsGroup
[
_fieldName_
](
i
,
dof
)
=
fieldsAtIPs
[
_fieldName_
](
*
iter
,
dof
);
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
gradFieldsGroup
[
_fieldName_
][
j
](
i
,
dof
)
=
gradFieldsAtIPs
[
_fieldName_
][
j
](
*
iter
,
dof
);
}
}
}
}
// calculate forces, & assemble
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
_fieldItr_
=
fields
.
find
(
_fieldName_
);
int
index
=
(
int
)
_fieldName_
;
if
(
physicsModel
->
null
(
_fieldName_
,
imat
))
continue
;
if
(
rhsMask
(
index
,
SOURCE
)
)
{
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
ndof
=
field
.
nCols
();
bool
has
=
physicsModel
->
weak_equation
(
_fieldName_
)
->
N_integrand
(
fieldsGroup
,
gradFieldsGroup
,
mat
,
group_Nfluxes
);
if
(
!
has
)
throw
ATC_Error
(
"atomic source can not be null currently"
);
int
i
=
0
;
for
(
set
<
int
>::
const_iterator
iter
=
indices
.
begin
();
iter
!=
indices
.
end
();
iter
++
,
i
++
)
{
for
(
int
dof
=
0
;
dof
<
ndof
;
++
dof
)
{
Nfluxes
[
_fieldName_
](
*
iter
,
dof
)
+=
group_Nfluxes
(
i
,
dof
);
}
}
}
}
// calculate forces, & assemble
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
physicsModel
->
null
(
_fieldName_
,
imat
))
continue
;
int
index
=
(
int
)
_fieldName_
;
if
(
rhsMask
(
index
,
FLUX
)
)
{
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
ndof
=
field
.
nCols
();
physicsModel
->
weak_equation
(
_fieldName_
)
->
B_integrand
(
fieldsGroup
,
gradFieldsGroup
,
mat
,
group_Bfluxes
);
int
i
=
0
;
for
(
set
<
int
>::
const_iterator
iter
=
indices
.
begin
();
iter
!=
indices
.
end
();
iter
++
,
i
++
)
{
for
(
int
dof
=
0
;
dof
<
ndof
;
++
dof
)
{
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
Bfluxes
[
_fieldName_
][
j
](
*
iter
,
dof
)
+=
group_Bfluxes
[
j
](
i
,
dof
);
}
}
}
}
}
}
}
// endif multiple materials
// assemble : N/Bfluxes -> rhs_local
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
int
index
=
(
int
)
_fieldName_
;
if
(
rhsMask
(
index
,
SOURCE
)
&&
Nfluxes
[
_fieldName_
].
nCols
()
>
0
)
{
rhs_local
[
_fieldName_
]
+=
N
.
transMat
(
weights
*
Nfluxes
[
_fieldName_
]);
}
if
(
rhsMask
(
index
,
FLUX
)
&&
(
Bfluxes
[
_fieldName_
][
0
]).
nCols
()
>
0
)
{
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
rhs_local
[
_fieldName_
]
+=
dN
[
j
]
->
transMat
(
weights
*
Bfluxes
[
_fieldName_
][
j
]);
}
}
}
}
// end nips > 0
// Share information between processors: rhs_local -> rhs
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
rhsMask
(
_fieldName_
,
FLUX
)
||
rhsMask
(
_fieldName_
,
SOURCE
))
{
DENS_MAT
&
myRhs
(
rhs
[
_fieldName_
].
set_quantity
());
int
count
=
rhs_local
[
_fieldName_
].
size
();
allsum
(
communicator_
,
rhs_local
[
_fieldName_
].
ptr
(),
myRhs
.
ptr
(),
count
);
}
}
}
// -------------------------------------------------------------
// compute sources using given quadrature
// -------------------------------------------------------------
void
FE_Engine
::
compute_source
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
set
<
int
>
>&
pointMaterialGroups
,
const
DIAG_MAT
&
weights
,
const
SPAR_MAT
&
N
,
const
SPAR_MAT_VEC
&
dN
,
FIELD_MATS
&
sources
)
const
{
int
nips
=
weights
.
nCols
();
if
(
nips
>
0
)
{
FIELD_MATS
Nfluxes
;
// compute fields and gradients of fields at given ips
GRAD_FIELD_MATS
gradFieldsAtIPs
;
FIELD_MATS
fieldsAtIPs
;
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
_fieldName_
=
_fieldItr_
->
first
;
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
dof
=
field
.
nCols
();
gradFieldsAtIPs
[
_fieldName_
].
assign
(
nSD_
,
DENS_MAT
(
nips
,
dof
));
fieldsAtIPs
[
_fieldName_
]
=
N
*
field
;
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
gradFieldsAtIPs
[
_fieldName_
][
j
]
=
(
*
dN
[
j
])
*
field
;
}
}
// treat single material point sets specially
int
nMatls
=
pointMaterialGroups
.
size
();
int
atomMatls
=
0
;
for
(
int
imat
=
0
;
imat
<
nMatls
;
imat
++
)
{
const
set
<
int
>
&
indices
=
pointMaterialGroups
(
imat
);
if
(
indices
.
size
()
>
0
)
atomMatls
++
;
}
bool
singleMaterial
=
(
atomMatls
==
1
);
// evaluate physics model
if
(
singleMaterial
)
{
int
imat
=
0
;
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
int
index
=
(
int
)
_fieldName_
;
if
(
!
physicsModel
->
null
(
_fieldName_
,
imat
))
{
if
(
rhsMask
(
index
,
SOURCE
)
)
{
bool
has
=
physicsModel
->
weak_equation
(
_fieldName_
)
->
N_integrand
(
fieldsAtIPs
,
gradFieldsAtIPs
,
mat
,
Nfluxes
[
_fieldName_
]);
if
(
!
has
)
throw
ATC_Error
(
"atomic source can not be null currently"
);
}
}
}
}
else
{
throw
ATC_Error
(
"compute source does not handle multiple materials currently"
);
}
// assemble
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
int
index
=
(
int
)
_fieldName_
;
if
(
rhsMask
(
index
,
SOURCE
)
)
{
sources
[
_fieldName_
]
=
weights
*
Nfluxes
[
_fieldName_
];
}
}
}
// no need to share information between processors
}
// -------------------------------------------------------------
// compute flux for post processing
// -------------------------------------------------------------
void
FE_Engine
::
compute_flux
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
int
>
&
elementMaterials
,
GRAD_FIELD_MATS
&
fluxes
,
const
DenseMatrix
<
bool
>
*
elementMask
)
const
{
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
rhsMask
(
_fieldName_
,
FLUX
))
{
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
fluxes
[
_fieldName_
].
assign
(
nSD_
,
DENS_MAT
(
field
.
nRows
(),
field
.
nCols
()));
}
}
// element wise assembly
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
// if element is masked, skip it
if
(
elementMask
&&
!
(
*
elementMask
)(
ielem
,
0
))
continue
;
// material id
int
imat
=
elementMaterials
(
ielem
);
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
interpolate_fields
(
ielem
,
fields
,
_conn_
,
_N_
,
_dN_
,
_weights_
,
_fieldsAtIPs_
,
_gradFieldsAtIPs_
);
_Nw_
=
_weights_
*
_N_
;
// evaluate Physics model & assemble
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
!
rhsMask
(
_fieldName_
,
FLUX
))
continue
;
if
(
physicsModel
->
null
(
_fieldName_
,
imat
))
continue
;
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
dof
=
field
.
nCols
();
physicsModel
->
weak_equation
(
_fieldName_
)
->
B_integrand
(
_fieldsAtIPs_
,
_gradFieldsAtIPs_
,
mat
,
_Bfluxes_
);
for
(
int
j
=
0
;
j
<
nSD_
;
j
++
)
{
_Nmat_
=
_Nw_
.
transMat
(
_Bfluxes_
[
j
]);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
k
=
0
;
k
<
dof
;
++
k
)
{
fluxes
[
_fieldName_
][
j
](
inode
,
k
)
+=
_Nmat_
(
i
,
k
);
}
}
}
}
}
// Assemble partial results
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
!
rhsMask
(
_fieldName_
,
FLUX
))
continue
;
for
(
int
j
=
0
;
j
<
nSD_
;
j
++
)
{
allsum
(
communicator_
,
MPI_IN_PLACE
,
fluxes
[
_fieldName_
][
j
].
ptr
(),
fluxes
[
_fieldName_
][
j
].
size
());
}
}
}
//-----------------------------------------------------------------
// boundary flux using native quadrature
//-----------------------------------------------------------------
void
FE_Engine
::
compute_boundary_flux
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
int
>
&
elementMaterials
,
const
set
<
pair
<
int
,
int
>
>
&
faceSet
,
FIELDS
&
rhs
)
const
{
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
rhsMask
(
_fieldName_
,
FLUX
))
{
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
rhs
[
_fieldName_
].
reset
(
field
.
nRows
(),
field
.
nCols
());
}
}
FIELD_MATS
localElementFields
;
GRAD_FIELD_MATS
gradFieldsAtIPs
;
FIELD_MATS
fieldsAtIPs
;
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
_fieldName_
=
_fieldItr_
->
first
;
const
FIELD
&
field
=
_fieldItr_
->
second
;
int
dof
=
field
.
nCols
();
gradFieldsAtIPs
[
_fieldName_
].
reserve
(
nSD_
);
for
(
int
i
=
0
;
i
<
nSD_
;
++
i
)
{
gradFieldsAtIPs
[
_fieldName_
].
push_back
(
DENS_MAT
(
nIPsPerFace_
,
dof
));
}
fieldsAtIPs
[
_fieldName_
].
reset
(
nIPsPerFace_
,
dof
);
localElementFields
[
_fieldName_
].
reset
(
nNodesPerElement_
,
dof
);
}
// element wise assembly
set
<
PAIR
>::
iterator
iter
;
for
(
iter
=
faceSet
.
begin
();
iter
!=
faceSet
.
end
();
iter
++
)
{
// get connectivity
int
ielem
=
iter
->
first
;
// if this is not our element, do not do calculations
if
(
!
feMesh_
->
is_owned_elt
(
ielem
))
continue
;
int
imat
=
elementMaterials
(
ielem
);
const
Material
*
mat
=
physicsModel
->
material
(
imat
);
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
// evaluate shape functions at ips
feMesh_
->
face_shape_function
(
*
iter
,
_fN_
,
_fdN_
,
_nN_
,
_fweights_
);
// interpolate fields and gradients of fields ips of this element
for
(
_fieldItr_
=
fields
.
begin
();
_fieldItr_
!=
fields
.
end
();
_fieldItr_
++
)
{
_fieldName_
=
_fieldItr_
->
first
;
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
dof
=
field
.
nCols
();
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
i
++
)
{
for
(
int
j
=
0
;
j
<
dof
;
j
++
)
{
localElementFields
[
_fieldName_
](
i
,
j
)
=
field
(
_conn_
(
i
),
j
);
}
}
// ips X dof = ips X nodes * nodes X dof
fieldsAtIPs
[
_fieldName_
]
=
_fN_
*
localElementFields
[
_fieldName_
];
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
gradFieldsAtIPs
[
_fieldName_
][
j
]
=
_fdN_
[
j
]
*
localElementFields
[
_fieldName_
];
}
}
// Evaluate- physics model
// do nothing for N_integrand
// nN : precomputed and held by ATC_Transfer
// assemble
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
int
index
=
(
int
)
_fieldName_
;
if
(
rhsMask
(
index
,
FLUX
)
)
{
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
DENS_MAT
&
myRhs
(
rhs
[
_fieldName_
].
set_quantity
());
physicsModel
->
weak_equation
(
_fieldName_
)
->
B_integrand
(
fieldsAtIPs
,
gradFieldsAtIPs
,
mat
,
_Bfluxes_
);
int
dof
=
field
.
nCols
();
for
(
int
j
=
0
;
j
<
nSD_
;
j
++
)
{
_Nmat_
=
_nN_
[
j
].
transMat
(
_fweights_
*
_Bfluxes_
[
j
]);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
k
=
0
;
k
<
dof
;
++
k
)
{
myRhs
(
inode
,
k
)
+=
_Nmat_
(
i
,
k
);
}
}
}
}
}
}
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
int
index
=
(
int
)
_fieldName_
;
if
(
rhsMask
(
index
,
FLUX
)
)
{
DENS_MAT
&
myRhs
(
rhs
[
_fieldName_
].
set_quantity
());
allsum
(
communicator_
,
MPI_IN_PLACE
,
myRhs
.
ptr
(),
myRhs
.
size
());
}
}
}
// -------------------------------------------------------------
// compute boundary flux using given quadrature and interpolation
// -------------------------------------------------------------
void
FE_Engine
::
compute_boundary_flux
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
PhysicsModel
*
physicsModel
,
const
Array
<
int
>
&
elementMaterials
,
const
Array
<
set
<
int
>
>
&
pointMaterialGroups
,
const
DIAG_MAT
&
_weights_
,
const
SPAR_MAT
&
N
,
const
SPAR_MAT_VEC
&
dN
,
const
DIAG_MAT
&
flux_mask
,
FIELDS
&
flux
,
const
DenseMatrix
<
bool
>
*
elementMask
,
const
set
<
int
>
*
nodeSet
)
const
{
FIELDS
rhs
,
rhsSubset
;
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
rhsMask
(
_fieldName_
,
FLUX
))
{
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
nrows
=
field
.
nRows
();
int
ncols
=
field
.
nCols
();
rhs
[
_fieldName_
].
reset
(
nrows
,
ncols
);
rhsSubset
[
_fieldName_
].
reset
(
nrows
,
ncols
);
}
}
// R_I = - int_Omega Delta _N_I .q dV
compute_rhs_vector
(
rhsMask
,
fields
,
physicsModel
,
elementMaterials
,
rhs
,
elementMask
);
// R_I^md = - int_Omega^md Delta _N_I .q dV
compute_rhs_vector
(
rhsMask
,
fields
,
physicsModel
,
pointMaterialGroups
,
_weights_
,
N
,
dN
,
rhsSubset
);
// flux_I = 1/Delta V_I V_I^md R_I + R_I^md
// where : Delta V_I = int_Omega _N_I dV
for
(
int
n
=
0
;
n
<
rhsMask
.
nRows
();
n
++
)
{
_fieldName_
=
FieldName
(
n
);
if
(
rhsMask
(
_fieldName_
,
FLUX
)
)
{
if
(
nodeSet
)
{
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
int
nrows
=
field
.
nRows
();
int
ncols
=
field
.
nCols
();
DENS_MAT
&
myFlux
(
flux
[
_fieldName_
].
set_quantity
());
const
DENS_MAT
&
myRhsSubset
(
rhsSubset
[
_fieldName_
].
quantity
());
const
DENS_MAT
&
myRhs
(
rhs
[
_fieldName_
].
quantity
());
myFlux
.
reset
(
nrows
,
ncols
);
set
<
int
>::
const_iterator
iset
;
for
(
iset
=
nodeSet
->
begin
();
iset
!=
nodeSet
->
end
();
iset
++
)
{
for
(
int
j
=
0
;
j
<
ncols
;
j
++
)
{
myFlux
(
*
iset
,
j
)
=
myRhsSubset
(
*
iset
,
j
)
-
flux_mask
(
*
iset
,
*
iset
)
*
myRhs
(
*
iset
,
j
);
}
}
}
else
{
flux
[
_fieldName_
]
=
rhsSubset
[
_fieldName_
].
quantity
()
-
flux_mask
*
(
rhs
[
_fieldName_
].
quantity
());
}
}
}
}
/** integrate a nodal field over an element set */
DENS_VEC
FE_Engine
::
integrate
(
const
DENS_MAT
&
field
,
const
ESET
&
eset
)
const
{
int
dof
=
field
.
nCols
();
DENS_MAT
eField
(
nNodesPerElement_
,
dof
);
int
nips
=
nIPsPerElement_
;
DENS_MAT
ipField
(
nips
,
dof
);
DENS_VEC
integral
(
dof
);
integral
=
0
;
for
(
ESET
::
const_iterator
itr
=
eset
.
begin
();
itr
!=
eset
.
end
();
++
itr
)
{
int
ielem
=
*
itr
;
// if this is not our element, do not do calculations
if
(
!
feMesh_
->
is_owned_elt
(
ielem
))
continue
;
feMesh_
->
shape_function
(
ielem
,
_N_
,
_weights_
);
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
i
++
)
{
for
(
int
j
=
0
;
j
<
dof
;
j
++
)
{
eField
(
i
,
j
)
=
field
(
_conn_
(
i
),
j
);
}}
ipField
=
_N_
*
eField
;
for
(
int
i
=
0
;
i
<
nips
;
++
i
)
{
for
(
int
j
=
0
;
j
<
dof
;
++
j
)
{
integral
(
j
)
+=
ipField
(
i
,
j
)
*
_weights_
[
i
];
}
}
}
// assemble partial results
allsum
(
communicator_
,
MPI_IN_PLACE
,
integral
.
ptr
(),
integral
.
size
());
return
integral
;
}
//-----------------------------------------------------------------
// Robin boundary flux using native quadrature
// integrate \int_delV _N_I s(x,t).n dA
//-----------------------------------------------------------------
void
FE_Engine
::
add_robin_fluxes
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
double
time
,
const
ROBIN_SURFACE_SOURCE
&
sourceFunctions
,
FIELDS
&
nodalSources
)
const
{
// sizing working arrays
DENS_MAT
xCoords
(
nSD_
,
nNodesPerElement_
);
DENS_MAT
faceSource
;
DENS_MAT
localField
;
DENS_MAT
xAtIPs
(
nSD_
,
nIPsPerFace_
);
DENS_MAT
uAtIPs
(
nIPsPerFace_
,
1
);
FIELDS
myNodalSources
;
// element wise assembly
ROBIN_SURFACE_SOURCE
::
const_iterator
src_iter
;
for
(
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
if
(
!
rhsMask
((
int
)
_fieldName_
,
ROBIN_SOURCE
))
continue
;
DENS_MAT
&
s
(
nodalSources
[
_fieldName_
].
set_quantity
());
myNodalSources
[
_fieldName_
]
=
DENS_MAT
(
s
.
nRows
(),
s
.
nCols
());
}
for
(
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
if
(
!
rhsMask
((
int
)
_fieldName_
,
ROBIN_SOURCE
))
continue
;
typedef
map
<
PAIR
,
Array
<
UXT_Function
*>
>
FSET
;
const
FSET
*
fset
=
(
const
FSET
*
)
&
(
src_iter
->
second
);
FSET
::
const_iterator
fset_iter
;
for
(
fset_iter
=
fset
->
begin
();
fset_iter
!=
fset
->
end
();
fset_iter
++
)
{
const
PAIR
&
face
=
fset_iter
->
first
;
const
int
elem
=
face
.
first
;
// if this is not our element, do not do calculations
if
(
!
feMesh_
->
is_owned_elt
(
elem
))
continue
;
const
Array
<
UXT_Function
*>
&
fs
=
fset_iter
->
second
;
_conn_
=
feMesh_
->
element_connectivity_unique
(
elem
);
// evaluate location at ips
feMesh_
->
face_shape_function
(
face
,
_fN_
,
_fdN_
,
_nN_
,
_fweights_
);
feMesh_
->
element_coordinates
(
elem
,
xCoords
);
xAtIPs
=
xCoords
*
(
_fN_
.
transpose
());
// collect field
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
feMesh_
->
element_field
(
elem
,
field
,
localField
);
uAtIPs
=
_fN_
*
localField
;
// interpolate prescribed flux at ips of this element
FSET
::
const_iterator
face_iter
=
fset
->
find
(
face
);
int
nFieldDOF
=
(
face_iter
->
second
).
size
();
faceSource
.
reset
(
nIPsPerFace_
,
nFieldDOF
);
for
(
int
ip
=
0
;
ip
<
nIPsPerFace_
;
++
ip
)
{
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
UXT_Function
*
f
=
fs
(
idof
);
if
(
!
f
)
continue
;
faceSource
(
ip
,
idof
)
=
f
->
f
(
&
(
uAtIPs
(
ip
,
0
)),
column
(
xAtIPs
,
ip
).
ptr
(),
time
);
DENS_MAT
coefsAtIPs
(
nIPsPerFace_
,
1
);
coefsAtIPs
(
ip
,
idof
)
=
f
->
dfdu
(
&
(
uAtIPs
(
ip
,
0
)),
column
(
xAtIPs
,
ip
).
ptr
(),
time
);
faceSource
(
ip
,
idof
)
-=
coefsAtIPs
(
ip
,
0
)
*
uAtIPs
(
ip
,
0
);
}
}
// assemble
DENS_MAT
&
s
(
myNodalSources
[
_fieldName_
].
set_quantity
());
_Nmat_
=
_fN_
.
transMat
(
_fweights_
*
faceSource
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
s
(
inode
,
idof
)
+=
_Nmat_
(
i
,
idof
);
}
}
}
}
// assemble partial result matrices
for
(
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
if
(
!
rhsMask
((
int
)
_fieldName_
,
ROBIN_SOURCE
))
continue
;
DENS_MAT
&
s
(
myNodalSources
[
_fieldName_
].
set_quantity
());
allsum
(
communicator_
,
MPI_IN_PLACE
,
s
.
ptr
(),
s
.
size
());
DENS_MAT
&
src
(
nodalSources
[
_fieldName_
].
set_quantity
());
src
+=
s
;
}
}
//-----------------------------------------------------------------
// Robin boundary flux stiffness using native quadrature
// integrate \int_delV _N_I ds/du(x,t).n dA
//-----------------------------------------------------------------
void
FE_Engine
::
add_robin_tangent
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
double
time
,
const
ROBIN_SURFACE_SOURCE
&
sourceFunctions
,
SPAR_MAT
&
tangent
)
const
{
// sizing working arrays
DENS_MAT
xCoords
(
nSD_
,
nNodesPerElement_
);
DENS_MAT
coefsAtIPs
;
DENS_MAT
localField
;
DENS_MAT
xAtIPs
(
nSD_
,
nIPsPerFace_
);
DENS_MAT
uAtIPs
(
nIPsPerFace_
,
1
);
// element wise assembly
ROBIN_SURFACE_SOURCE
::
const_iterator
src_iter
;
SPAR_MAT
K
(
tangent
.
nRows
(),
tangent
.
nCols
());
for
(
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
if
(
!
rhsMask
((
int
)
_fieldName_
,
ROBIN_SOURCE
))
continue
;
typedef
map
<
PAIR
,
Array
<
UXT_Function
*>
>
FSET
;
const
FSET
*
fset
=
(
const
FSET
*
)
&
(
src_iter
->
second
);
FSET
::
const_iterator
fset_iter
;
for
(
fset_iter
=
fset
->
begin
();
fset_iter
!=
fset
->
end
();
fset_iter
++
)
{
const
PAIR
&
face
=
fset_iter
->
first
;
const
int
elem
=
face
.
first
;
// if this is not our element, do not do calculations
if
(
!
feMesh_
->
is_owned_elt
(
elem
))
continue
;
const
Array
<
UXT_Function
*>
&
fs
=
fset_iter
->
second
;
_conn_
=
feMesh_
->
element_connectivity_unique
(
elem
);
// evaluate location at ips
feMesh_
->
face_shape_function
(
face
,
_fN_
,
_fdN_
,
_nN_
,
_fweights_
);
feMesh_
->
element_coordinates
(
elem
,
xCoords
);
xAtIPs
=
xCoords
*
(
_fN_
.
transpose
());
// collect field
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
feMesh_
->
element_field
(
elem
,
field
,
localField
);
uAtIPs
=
_fN_
*
localField
;
// interpolate prescribed flux at ips of this element
FSET
::
const_iterator
face_iter
=
fset
->
find
(
face
);
int
nFieldDOF
=
(
face_iter
->
second
).
size
();
coefsAtIPs
.
reset
(
nIPsPerFace_
,
nFieldDOF
);
for
(
int
ip
=
0
;
ip
<
nIPsPerFace_
;
++
ip
)
{
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
UXT_Function
*
f
=
fs
(
idof
);
if
(
!
f
)
continue
;
coefsAtIPs
(
ip
,
idof
)
=
f
->
dfdu
(
&
(
uAtIPs
(
ip
,
0
)),
column
(
xAtIPs
,
ip
).
ptr
(),
time
);
}
}
// assemble
DIAG_MAT
D
(
column
(
coefsAtIPs
,
0
));
D
*=
-
1
;
D
*=
_fweights_
;
_Nmat_
=
_fN_
.
transMat
(
D
*
_fN_
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
_conn_
(
j
);
K
.
add
(
inode
,
jnode
,
_Nmat_
(
i
,
j
));
}
}
}
}
// assemble partial result matrices
#ifdef ISOLATE_FE
sparse_allsum
(
communicator_
,
K
);
#else
LammpsInterface
::
instance
()
->
sparse_allsum
(
K
);
#endif
tangent
+=
K
;
}
//-----------------------------------------------------------------
// Robin boundary flux using native quadrature
// integrate \int_delV _N_I s(x,t).n dA
//-----------------------------------------------------------------
void
FE_Engine
::
add_open_fluxes
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
OPEN_SURFACE
&
openFaces
,
FIELDS
&
nodalSources
,
const
FieldName
Velocity
)
const
{
// sizing working arrays
DENS_MAT
xCoords
(
nSD_
,
nNodesPerElement_
);
DENS_MAT
faceSource
;
DENS_MAT
localField
;
DENS_MAT
xAtIPs
(
nSD_
,
nIPsPerFace_
);
DENS_MAT
uAtIPs
;
// (nIPsPerFace_,field nCols);
DENS_MAT
aAtIPs
(
nIPsPerFace_
,
1
);
FIELDS
myNodalSources
;
// throw error if electron velocity is not defined
// element wise assembly
OPEN_SURFACE
::
const_iterator
face_iter
;
for
(
face_iter
=
openFaces
.
begin
();
face_iter
!=
openFaces
.
end
();
face_iter
++
)
{
_fieldName_
=
face_iter
->
first
;
if
(
!
rhsMask
((
int
)
_fieldName_
,
OPEN_SOURCE
))
continue
;
DENS_MAT
&
s
(
nodalSources
[
_fieldName_
].
set_quantity
());
myNodalSources
[
_fieldName_
]
=
DENS_MAT
(
s
.
nRows
(),
s
.
nCols
());
}
for
(
face_iter
=
openFaces
.
begin
();
face_iter
!=
openFaces
.
end
();
face_iter
++
)
{
_fieldName_
=
face_iter
->
first
;
if
(
!
rhsMask
((
int
)
_fieldName_
,
OPEN_SOURCE
))
continue
;
typedef
set
<
PAIR
>
FSET
;
const
FSET
*
fset
=
(
const
FSET
*
)
&
(
face_iter
->
second
);
FSET
::
const_iterator
fset_iter
;
for
(
fset_iter
=
fset
->
begin
();
fset_iter
!=
fset
->
end
();
fset_iter
++
)
{
const
PAIR
&
face
=
*
fset_iter
;
const
int
elem
=
face
.
first
;
// if this is not our element, do not do calculations
if
(
!
feMesh_
->
is_owned_elt
(
elem
))
continue
;
// get velocity, multiply with field (vector gives tensor)
// dot with face normal
_conn_
=
feMesh_
->
element_connectivity_unique
(
elem
);
// evaluate location at ips
feMesh_
->
face_shape_function
(
face
,
_fN_
,
_fdN_
,
_nN_
,
_fweights_
);
// collect field at ips
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
feMesh_
->
element_field
(
elem
,
field
,
localField
);
int
nFieldDOF
=
field
.
nCols
();
// "u" is the quantity being advected
uAtIPs
.
reset
(
nIPsPerFace_
,
nFieldDOF
);
uAtIPs
=
_fN_
*
localField
;
// collect velocity at ips for "advection" = v.n
_fieldItr_
=
fields
.
find
(
Velocity
);
const
DENS_MAT
&
advection
=
(
_fieldItr_
->
second
).
quantity
();
// advection velocity
feMesh_
->
element_field
(
elem
,
advection
,
localField
);
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
// nSD_ == nDOF for the velocity
for
(
int
ip
=
0
;
ip
<
nIPsPerFace_
;
++
ip
)
{
for
(
int
I
=
0
;
I
<
nNodesPerElement_
;
++
I
)
{
aAtIPs
(
ip
,
0
)
+=
(
_nN_
[
j
])(
ip
,
I
)
*
localField
(
I
,
j
);
}
}
}
// construct open flux at ips of this element
// flux = field \otimes advection_vector \dot n
faceSource
.
reset
(
nIPsPerFace_
,
nFieldDOF
);
for
(
int
ip
=
0
;
ip
<
nIPsPerFace_
;
++
ip
)
{
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
// multiply field DOF value times advection vector
faceSource
(
ip
,
idof
)
=
aAtIPs
(
ip
,
1
)
*
uAtIPs
(
ip
,
idof
);
//(v.n) u
}
}
// assemble contributions for this direction of the face normal
// _nN_[j](ip,I) nodal shapefunction(I) at ip X face normal(j)
// _fweights_ diagonal nips X nips ==> facesource is nips X ndofs
DENS_MAT
&
s
(
myNodalSources
[
_fieldName_
].
set_quantity
());
_Nmat_
=
_fN_
.
transMat
(
_fweights_
*
faceSource
);
// s_Ii = \sum_ip N_I (u_i v.n)_ip wg_ip
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
s
(
inode
,
idof
)
+=
_Nmat_
(
i
,
idof
);
}
}
}
}
// assemble partial result matrices
for
(
face_iter
=
openFaces
.
begin
();
face_iter
!=
openFaces
.
end
();
face_iter
++
)
{
_fieldName_
=
face_iter
->
first
;
if
(
!
rhsMask
((
int
)
_fieldName_
,
OPEN_SOURCE
))
continue
;
DENS_MAT
&
s
(
myNodalSources
[
_fieldName_
].
set_quantity
());
allsum
(
communicator_
,
MPI_IN_PLACE
,
s
.
ptr
(),
s
.
size
());
DENS_MAT
&
src
(
nodalSources
[
_fieldName_
].
set_quantity
());
src
+=
s
;
}
}
//-----------------------------------------------------------------
// Open boundary flux stiffness using native quadrature
// integrate \int_delV _N_I ds/du(x,t).n dA
//-----------------------------------------------------------------
void
FE_Engine
::
add_open_tangent
(
const
Array2D
<
bool
>
&
rhsMask
,
const
FIELDS
&
fields
,
const
OPEN_SURFACE
&
openFaces
,
SPAR_MAT
&
tangent
,
const
FieldName
Velocity
)
const
{
// sizing working arrays
DENS_MAT
xCoords
(
nSD_
,
nNodesPerElement_
);
DENS_MAT
faceSource
;
DENS_MAT
localField
;
DENS_MAT
xAtIPs
(
nSD_
,
nIPsPerFace_
);
DENS_MAT
uAtIPs
;
// (nIPsPerFace_,field nCols);
DENS_MAT
aAtIPs
(
nIPsPerFace_
,
nSD_
);
SPAR_MAT
K
(
tangent
.
nRows
(),
tangent
.
nCols
());
// element wise assembly
OPEN_SURFACE
::
const_iterator
face_iter
;
for
(
face_iter
=
openFaces
.
begin
();
face_iter
!=
openFaces
.
end
();
face_iter
++
)
{
_fieldName_
=
face_iter
->
first
;
if
(
!
rhsMask
((
int
)
_fieldName_
,
OPEN_SOURCE
))
continue
;
bool
convective
=
false
;
if
(
_fieldName_
==
Velocity
)
convective
=
true
;
typedef
set
<
PAIR
>
FSET
;
const
FSET
*
fset
=
(
const
FSET
*
)
&
(
face_iter
->
second
);
FSET
::
const_iterator
fset_iter
;
for
(
fset_iter
=
fset
->
begin
();
fset_iter
!=
fset
->
end
();
fset_iter
++
)
{
const
PAIR
&
face
=
*
fset_iter
;
const
int
elem
=
face
.
first
;
// if this is not our element, do not do calculations
if
(
!
feMesh_
->
is_owned_elt
(
elem
))
continue
;
_conn_
=
feMesh_
->
element_connectivity_unique
(
elem
);
// evaluate location at ips
feMesh_
->
face_shape_function
(
face
,
_fN_
,
_fdN_
,
_nN_
,
_fweights_
);
// collect field at ips
_fieldItr_
=
fields
.
find
(
_fieldName_
);
const
DENS_MAT
&
field
=
(
_fieldItr_
->
second
).
quantity
();
feMesh_
->
element_field
(
elem
,
field
,
localField
);
int
nFieldDOF
=
field
.
nCols
();
// "u" is the quantity being advected
uAtIPs
.
reset
(
nIPsPerFace_
,
nFieldDOF
);
uAtIPs
=
_fN_
*
localField
;
// collect velocity at ips for "advection" = v.n
_fieldItr_
=
fields
.
find
(
Velocity
);
const
DENS_MAT
&
advection
=
(
_fieldItr_
->
second
).
quantity
();
// advection velocity
feMesh_
->
element_field
(
elem
,
advection
,
localField
);
for
(
int
j
=
0
;
j
<
nSD_
;
++
j
)
{
// nSD_ == nDOF for the velocity
for
(
int
ip
=
0
;
ip
<
nIPsPerFace_
;
++
ip
)
{
for
(
int
I
=
0
;
I
<
nNodesPerElement_
;
++
I
)
{
aAtIPs
(
ip
,
0
)
+=
(
_nN_
[
j
])(
ip
,
I
)
*
localField
(
I
,
j
);
}
}
}
// K_IJ = \sum_ip N_I ( v.n I + u (x) n ) N_J wg_ip
DENS_MAT
D
(
nFieldDOF
,
nFieldDOF
);
DENS_VEC
n
(
nSD_
);
for
(
int
ip
=
0
;
ip
<
nIPsPerFace_
;
++
ip
)
{
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
D
(
idof
,
idof
)
-=
aAtIPs
(
ip
,
0
);
if
(
convective
)
{
feMesh_
->
face_normal
(
face
,
ip
,
n
);
for
(
int
jdof
=
0
;
jdof
<
nFieldDOF
;
++
jdof
)
{
D
(
idof
,
jdof
)
-=
uAtIPs
(
ip
,
idof
)
*
n
(
jdof
);
}
}
}
}
// assemble
_Nmat_
=
_fN_
.
transMat
(
D
*
_fN_
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
_conn_
(
j
);
K
.
add
(
inode
,
jnode
,
_Nmat_
(
i
,
j
));
}
}
}
}
// assemble partial result matrices
#ifdef ISOLATE_FE
sparse_allsum
(
communicator_
,
K
);
#else
LammpsInterface
::
instance
()
->
sparse_allsum
(
K
);
#endif
tangent
+=
K
;
}
//-----------------------------------------------------------------
// prescribed boundary flux using native quadrature
// integrate \int_delV _N_I s(x,t).n dA
//-----------------------------------------------------------------
void
FE_Engine
::
add_fluxes
(
const
Array
<
bool
>
&
fieldMask
,
const
double
time
,
const
SURFACE_SOURCE
&
sourceFunctions
,
FIELDS
&
nodalSources
)
const
{
// sizing working arrays
DENS_MAT
xCoords
(
nSD_
,
nNodesPerElement_
);
DENS_MAT
xAtIPs
(
nSD_
,
nIPsPerFace_
);
DENS_MAT
faceSource
;
FIELDS
myNodalSources
;
// element wise assembly
SURFACE_SOURCE
::
const_iterator
src_iter
;
for
(
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
if
(
!
fieldMask
((
int
)
_fieldName_
))
continue
;
DENS_MAT
&
s
(
nodalSources
[
_fieldName_
].
set_quantity
());
myNodalSources
[
_fieldName_
]
=
DENS_MAT
(
s
.
nRows
(),
s
.
nCols
());
}
for
(
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
if
(
!
fieldMask
((
int
)
_fieldName_
))
continue
;
typedef
map
<
PAIR
,
Array
<
XT_Function
*>
>
FSET
;
const
FSET
*
fset
=
(
const
FSET
*
)
&
(
src_iter
->
second
);
FSET
::
const_iterator
fset_iter
;
for
(
fset_iter
=
fset
->
begin
();
fset_iter
!=
fset
->
end
();
fset_iter
++
)
{
const
PAIR
&
face
=
fset_iter
->
first
;
const
int
elem
=
face
.
first
;
// if this is not our element, do not do calculations
if
(
!
feMesh_
->
is_owned_elt
(
elem
))
continue
;
const
Array
<
XT_Function
*>
&
fs
=
fset_iter
->
second
;
_conn_
=
feMesh_
->
element_connectivity_unique
(
elem
);
// evaluate location at ips
feMesh_
->
face_shape_function
(
face
,
_fN_
,
_fdN_
,
_nN_
,
_fweights_
);
feMesh_
->
element_coordinates
(
elem
,
xCoords
);
MultAB
(
xCoords
,
_fN_
,
xAtIPs
,
0
,
1
);
//xAtIPs = xCoords*(N.transpose());
// interpolate prescribed flux at ips of this element
FSET
::
const_iterator
face_iter
=
fset
->
find
(
face
);
if
(
face_iter
==
fset
->
end
())
{
stringstream
ss
;
ss
<<
"face not found"
<<
std
::
endl
;
print_msg
(
communicator_
,
ss
.
str
());
}
int
nFieldDOF
=
(
face_iter
->
second
).
size
();
faceSource
.
reset
(
nIPsPerFace_
,
nFieldDOF
);
for
(
int
ip
=
0
;
ip
<
nIPsPerFace_
;
++
ip
)
{
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
XT_Function
*
f
=
fs
(
idof
);
if
(
!
f
)
continue
;
faceSource
(
ip
,
idof
)
=
f
->
f
(
column
(
xAtIPs
,
ip
).
ptr
(),
time
);
}
}
// assemble
DENS_MAT
&
s
(
myNodalSources
[
_fieldName_
].
set_quantity
());
_Nmat_
=
_fN_
.
transMat
(
_fweights_
*
faceSource
);
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
s
(
inode
,
idof
)
+=
_Nmat_
(
i
,
idof
);
}
// end assemble nFieldDOF
}
// end assemble nNodesPerElement_
}
// end fset loop
}
// field loop
// assemble partial result matrices
for
(
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
if
(
!
fieldMask
((
int
)
_fieldName_
))
continue
;
DENS_MAT
&
s
(
myNodalSources
[
_fieldName_
].
set_quantity
());
allsum
(
communicator_
,
MPI_IN_PLACE
,
s
.
ptr
(),
s
.
size
());
DENS_MAT
&
src
(
nodalSources
[
_fieldName_
].
set_quantity
());
src
+=
s
;
}
}
//-----------------------------------------------------------------
// prescribed volume flux using native quadrature
// integrate \int_V _N_I s(x,t) dV
//-----------------------------------------------------------------
void
FE_Engine
::
add_sources
(
const
Array
<
bool
>
&
fieldMask
,
const
double
time
,
const
VOLUME_SOURCE
&
sourceFunctions
,
FIELDS
&
nodalSources
)
const
{
// sizing working arrays
DENS_MAT
elemSource
;
DENS_MAT
xCoords
(
nSD_
,
nNodesPerElement_
);
DENS_MAT
xAtIPs
(
nSD_
,
nIPsPerElement_
);
FIELDS
myNodalSources
;
for
(
VOLUME_SOURCE
::
const_iterator
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
int
index
=
(
int
)
_fieldName_
;
if
(
!
fieldMask
(
index
))
continue
;
DENS_MAT
&
s
(
nodalSources
[
_fieldName_
].
set_quantity
());
myNodalSources
[
_fieldName_
]
=
DENS_MAT
(
s
.
nRows
(),
s
.
nCols
());
}
vector
<
int
>
myElems
=
feMesh_
->
owned_elts
();
for
(
vector
<
int
>::
iterator
elemsIter
=
myElems
.
begin
();
elemsIter
!=
myElems
.
end
();
++
elemsIter
)
{
int
ielem
=
*
elemsIter
;
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
// evaluate location at ips
feMesh_
->
shape_function
(
ielem
,
_N_
,
_weights_
);
feMesh_
->
element_coordinates
(
ielem
,
xCoords
);
xAtIPs
=
xCoords
*
(
_N_
.
transpose
());
for
(
VOLUME_SOURCE
::
const_iterator
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
int
index
=
(
int
)
_fieldName_
;
if
(
fieldMask
(
index
)
)
{
const
Array2D
<
XT_Function
*>
*
thisSource
=
(
const
Array2D
<
XT_Function
*>
*
)
&
(
src_iter
->
second
);
int
nFieldDOF
=
thisSource
->
nCols
();
elemSource
.
reset
(
nIPsPerElement_
,
nFieldDOF
);
// interpolate prescribed flux at ips of this element
for
(
int
ip
=
0
;
ip
<
nIPsPerElement_
;
++
ip
)
{
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
XT_Function
*
f
=
(
*
thisSource
)(
ielem
,
idof
);
if
(
f
)
{
elemSource
(
ip
,
idof
)
=
f
->
f
(
column
(
xAtIPs
,
ip
).
ptr
(),
time
);
}
}
}
// assemble
_Nmat_
=
_N_
.
transMat
(
_weights_
*
elemSource
);
DENS_MAT
&
s
(
myNodalSources
[
_fieldName_
].
set_quantity
());
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
++
i
)
{
int
inode
=
_conn_
(
i
);
for
(
int
idof
=
0
;
idof
<
nFieldDOF
;
++
idof
)
{
s
(
inode
,
idof
)
+=
_Nmat_
(
i
,
idof
);
}
}
}
}
}
// Aggregate unmasked nodal sources on all processors.
for
(
VOLUME_SOURCE
::
const_iterator
src_iter
=
sourceFunctions
.
begin
();
src_iter
!=
sourceFunctions
.
end
();
src_iter
++
)
{
_fieldName_
=
src_iter
->
first
;
int
index
=
(
int
)
_fieldName_
;
if
(
!
fieldMask
(
index
))
continue
;
DENS_MAT
&
s
(
myNodalSources
[
_fieldName_
].
set_quantity
());
allsum
(
communicator_
,
MPI_IN_PLACE
,
s
.
ptr
(),
s
.
size
());
DENS_MAT
&
src
(
nodalSources
[
_fieldName_
].
set_quantity
());
src
+=
s
;
}
}
//-----------------------------------------------------------------
// previously computed nodal sources
//-----------------------------------------------------------------
void
FE_Engine
::
add_sources
(
const
Array
<
bool
>
&
fieldMask
,
const
double
time
,
const
FIELDS
&
sources
,
FIELDS
&
nodalSources
)
const
{
FIELDS
::
const_iterator
itr
;
for
(
itr
=
sources
.
begin
();
itr
!=
sources
.
end
();
itr
++
)
{
_fieldName_
=
itr
->
first
;
if
(
!
fieldMask
((
int
)
_fieldName_
))
continue
;
DENS_MAT
&
src
(
nodalSources
[
_fieldName_
].
set_quantity
());
const
DENS_MAT
&
s
((
sources
.
find
(
_fieldName_
)
->
second
).
quantity
());
for
(
int
inode
=
0
;
inode
<
nNodesUnique_
;
++
inode
)
{
for
(
int
j
=
0
;
j
<
src
.
nCols
();
++
j
)
{
src
(
inode
,
j
)
+=
s
(
inode
,
j
);
}
}
}
}
//-----------------------------------------------------------------
// boundary integral of a nodal field
//-----------------------------------------------------------------
void
FE_Engine
::
field_surface_flux
(
const
DENS_MAT
&
field
,
const
set
<
PAIR
>
&
faceSet
,
DENS_MAT
&
values
,
const
bool
contour
,
const
int
axis
)
const
{
int
dof
=
field
.
nCols
();
double
a
[
3
]
=
{
0
,
0
,
0
};
a
[
axis
]
=
1
;
// sizing working arrays
DENS_MAT
n
(
nSD_
,
nIPsPerFace_
);
DENS_MAT
localElementFields
(
nNodesPerElement_
,
dof
);
DENS_MAT
integrals
(
dof
,
nSD_
);
DENS_MAT
fieldsAtIPs
;
// SJL shouldn't this just be _fieldsAtIPs_
// the data member?
// element wise assembly
set
<
pair
<
int
,
int
>
>::
iterator
iter
;
for
(
iter
=
faceSet
.
begin
();
iter
!=
faceSet
.
end
();
iter
++
)
{
int
ielem
=
iter
->
first
;
// if this is not our element, do not do calculations
//if (!feMesh_->is_owned_elt(ielem)) continue;
// evaluate shape functions at ips
feMesh_
->
face_shape_function
(
*
iter
,
_N_
,
n
,
_fweights_
);
// cross n with axis to get tangent
if
(
contour
)
{
double
t
[
3
];
for
(
int
i
=
0
;
i
<
nIPsPerFace_
;
i
++
)
{
t
[
0
]
=
a
[
1
]
*
n
(
2
,
i
)
-
a
[
2
]
*
n
(
1
,
i
);
t
[
1
]
=
a
[
2
]
*
n
(
0
,
i
)
-
a
[
0
]
*
n
(
2
,
i
);
t
[
2
]
=
a
[
0
]
*
n
(
1
,
i
)
-
a
[
1
]
*
n
(
0
,
i
);
n
(
0
,
i
)
=
t
[
0
];
n
(
1
,
i
)
=
t
[
1
];
n
(
2
,
i
)
=
t
[
2
];
}
}
// get connectivity
_conn_
=
feMesh_
->
element_connectivity_unique
(
ielem
);
// interpolate fields and gradients of fields ips of this element
for
(
int
i
=
0
;
i
<
nNodesPerElement_
;
i
++
)
{
for
(
int
j
=
0
;
j
<
dof
;
j
++
)
{
localElementFields
(
i
,
j
)
=
field
(
_conn_
(
i
),
j
);
}
}
// ips X dof = ips X nodes * nodes X dof
fieldsAtIPs
=
_N_
*
localElementFields
;
// assemble : integral(k,j) = sum_ip n(j,ip) wg(ip,ip) field(ip,k)
_Nmat_
=
n
*
_fweights_
*
fieldsAtIPs
;
for
(
int
j
=
0
;
j
<
nSD_
;
j
++
)
{
for
(
int
k
=
0
;
k
<
dof
;
++
k
)
{
integrals
(
k
,
j
)
+=
_Nmat_
(
j
,
k
);
}
}
}
// assemble partial results
//LammpsInterface::instance()->allsum(MPI_IN_PLACE, integrals.ptr(), integrals.size());
// (S.n)_1 = S_1j n_j = S_11 n_1 + S_12 n_2 + S_13 n_3
// (S.n)_2 = S_2j n_j = S_21 n_1 + S_22 n_2 + S_23 n_3
// (S.n)_3 = S_3j n_j = S_31 n_1 + S_32 n_2 + S_33 n_3
if
(
dof
==
9
)
{
// tensor
values
.
reset
(
nSD_
,
1
);
values
(
0
,
0
)
=
integrals
(
0
,
0
)
+
integrals
(
1
,
1
)
+
integrals
(
2
,
2
);
values
(
1
,
0
)
=
integrals
(
3
,
0
)
+
integrals
(
4
,
1
)
+
integrals
(
5
,
2
);
values
(
2
,
0
)
=
integrals
(
6
,
0
)
+
integrals
(
7
,
1
)
+
integrals
(
8
,
2
);
}
else
if
(
dof
==
6
)
{
// sym tensor
values
.
reset
(
nSD_
,
1
);
values
(
0
,
0
)
=
integrals
(
0
,
0
)
+
integrals
(
3
,
1
)
+
integrals
(
4
,
2
);
values
(
1
,
0
)
=
integrals
(
3
,
0
)
+
integrals
(
1
,
1
)
+
integrals
(
5
,
2
);
values
(
2
,
0
)
=
integrals
(
4
,
1
)
+
integrals
(
5
,
1
)
+
integrals
(
2
,
2
);
}
// (v.n) = v_j n_j
else
if
(
dof
==
3
)
{
// vector
values
.
reset
(
1
,
1
);
values
(
0
,
0
)
=
integrals
(
0
,
0
)
+
integrals
(
1
,
1
)
+
integrals
(
2
,
2
);
}
// s n = s n_j e_j
else
if
(
dof
==
1
)
{
// scalar
values
.
reset
(
nSD_
,
1
);
values
(
0
,
0
)
=
integrals
(
0
,
0
);
values
(
1
,
0
)
=
integrals
(
0
,
1
);
values
(
2
,
0
)
=
integrals
(
0
,
2
);
}
else
{
string
msg
=
"FE_Engine::field_surface_flux unsupported field width: "
;
msg
+=
to_string
(
dof
);
throw
ATC_Error
(
msg
);
}
}
//-----------------------------------------------------------------
// evaluate shape functions at given points
//-----------------------------------------------------------------
void
FE_Engine
::
evaluate_shape_functions
(
const
MATRIX
&
pt_coords
,
SPAR_MAT
&
N
)
const
{
// Get shape function and derivatives at atomic locations
int
nnodes
=
feMesh_
->
num_nodes_unique
();
int
npts
=
pt_coords
.
nRows
();
// loop over point (atom) coordinates
DENS_VEC
x
(
nSD_
);
Array
<
int
>
node_index
(
nNodesPerElement_
);
DENS_VEC
shp
(
nNodesPerElement_
);
N
.
reset
(
npts
,
nnodes
);
for
(
int
i
=
0
;
i
<
npts
;
++
i
)
{
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
x
(
k
)
=
pt_coords
(
i
,
k
);
}
feMesh_
->
shape_functions
(
x
,
shp
,
node_index
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
node_index
(
j
);
N
.
add
(
i
,
jnode
,
shp
(
j
));
}
}
N
.
compress
();
}
//-----------------------------------------------------------------
// evaluate shape functions and their spatial derivatives at given points
//-----------------------------------------------------------------
void
FE_Engine
::
evaluate_shape_functions
(
const
MATRIX
&
pt_coords
,
SPAR_MAT
&
N
,
SPAR_MAT_VEC
&
dN
)
const
{
// Get shape function and derivatives at atomic locations
int
nnodes
=
feMesh_
->
num_nodes_unique
();
int
npts
=
pt_coords
.
nRows
();
// loop over point (atom) coordinates
DENS_VEC
x
(
nSD_
);
Array
<
int
>
node_index
(
nNodesPerElement_
);
DENS_VEC
shp
(
nNodesPerElement_
);
DENS_MAT
dshp
(
nSD_
,
nNodesPerElement_
);
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
dN
[
k
]
->
reset
(
npts
,
nnodes
);
}
N
.
reset
(
npts
,
nnodes
);
for
(
int
i
=
0
;
i
<
npts
;
++
i
)
{
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
x
(
k
)
=
pt_coords
(
i
,
k
);
}
feMesh_
->
shape_functions
(
x
,
shp
,
dshp
,
node_index
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
node_index
(
j
);
N
.
add
(
i
,
jnode
,
shp
(
j
));
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
dN
[
k
]
->
add
(
i
,
jnode
,
dshp
(
k
,
j
));
}
}
}
N
.
compress
();
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
dN
[
k
]
->
compress
();
}
}
//-----------------------------------------------------------------
// evaluate shape functions at given points
//-----------------------------------------------------------------
void
FE_Engine
::
evaluate_shape_functions
(
const
MATRIX
&
pt_coords
,
const
INT_ARRAY
&
pointToEltMap
,
SPAR_MAT
&
N
)
const
{
// Get shape function and derivatives at atomic locations
int
nnodes
=
feMesh_
->
num_nodes_unique
();
int
npts
=
pt_coords
.
nRows
();
// loop over point (atom) coordinates
DENS_VEC
x
(
nSD_
);
Array
<
int
>
node_index
(
nNodesPerElement_
);
DENS_VEC
shp
(
nNodesPerElement_
);
N
.
reset
(
npts
,
nnodes
);
for
(
int
i
=
0
;
i
<
npts
;
++
i
)
{
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
x
(
k
)
=
pt_coords
(
i
,
k
);
}
int
eltId
=
pointToEltMap
(
i
,
0
);
feMesh_
->
shape_functions
(
x
,
eltId
,
shp
,
node_index
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
node_index
(
j
);
N
.
add
(
i
,
jnode
,
shp
(
j
));
}
}
N
.
compress
();
}
//-----------------------------------------------------------------
// evaluate shape functions and their spatial derivatives at given points
//-----------------------------------------------------------------
void
FE_Engine
::
evaluate_shape_functions
(
const
MATRIX
&
pt_coords
,
const
INT_ARRAY
&
pointToEltMap
,
SPAR_MAT
&
N
,
SPAR_MAT_VEC
&
dN
)
const
{
// Get shape function and derivatives at atomic locations
int
nnodes
=
feMesh_
->
num_nodes_unique
();
int
npts
=
pt_coords
.
nRows
();
// loop over point (atom) coordinates
DENS_VEC
x
(
nSD_
);
Array
<
int
>
node_index
(
nNodesPerElement_
);
DENS_VEC
shp
(
nNodesPerElement_
);
DENS_MAT
dshp
(
nSD_
,
nNodesPerElement_
);
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
dN
[
k
]
->
reset
(
npts
,
nnodes
);
}
N
.
reset
(
npts
,
nnodes
);
for
(
int
i
=
0
;
i
<
npts
;
++
i
)
{
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
x
(
k
)
=
pt_coords
(
i
,
k
);
}
int
eltId
=
pointToEltMap
(
i
,
0
);
feMesh_
->
shape_functions
(
x
,
eltId
,
shp
,
dshp
,
node_index
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
node_index
(
j
);
N
.
add
(
i
,
jnode
,
shp
(
j
));
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
dN
[
k
]
->
add
(
i
,
jnode
,
dshp
(
k
,
j
));
}
}
}
N
.
compress
();
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
dN
[
k
]
->
compress
();
}
}
//-----------------------------------------------------------------
// evaluate shape functions spatial derivatives at given points
//-----------------------------------------------------------------
void
FE_Engine
::
evaluate_shape_function_derivatives
(
const
MATRIX
&
pt_coords
,
const
INT_ARRAY
&
pointToEltMap
,
SPAR_MAT_VEC
&
dNdx
)
const
{
int
nnodes
=
feMesh_
->
num_nodes_unique
();
int
npts
=
pt_coords
.
nRows
();
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
dNdx
[
k
]
->
reset
(
npts
,
nnodes
);
}
// loop over point (atom) coordinates
DENS_VEC
x
(
nSD_
);
Array
<
int
>
node_index
(
nNodesPerElement_
);
DENS_MAT
dshp
(
nSD_
,
nNodesPerElement_
);
for
(
int
i
=
0
;
i
<
npts
;
++
i
)
{
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
x
(
k
)
=
pt_coords
(
i
,
k
);
}
int
eltId
=
pointToEltMap
(
i
,
0
);
feMesh_
->
shape_function_derivatives
(
x
,
eltId
,
dshp
,
node_index
);
for
(
int
j
=
0
;
j
<
nNodesPerElement_
;
++
j
)
{
int
jnode
=
node_index
(
j
);
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
dNdx
[
k
]
->
add
(
i
,
jnode
,
dshp
(
k
,
j
));
}
}
}
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
dNdx
[
k
]
->
compress
();
}
}
//-------------------------------------------------------------------
// set kernels
//-------------------------------------------------------------------
void
FE_Engine
::
set_kernel
(
KernelFunction
*
ptr
){
kernelFunction_
=
ptr
;}
//-------------------------------------------------------------------
// kernel_matrix_bandwidth
//-------------------------------------------------------------------
int
FE_Engine
::
kernel_matrix_bandwidth
(
const
MATRIX
&
ptCoords
)
const
{
if
(
!
kernelFunction_
)
throw
ATC_Error
(
"no kernel function specified"
);
int
npts
=
ptCoords
.
nRows
();
int
bw
=
0
;
if
(
npts
>
0
)
{
DENS_VEC
xI
(
nSD_
),
x
(
nSD_
),
dx
(
nSD_
);
for
(
int
i
=
0
;
i
<
nNodes_
;
i
++
)
{
xI
=
feMesh_
->
global_coordinates
(
i
);
for
(
int
j
=
0
;
j
<
npts
;
j
++
)
{
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
x
(
k
)
=
ptCoords
(
j
,
k
);
}
dx
=
x
-
xI
;
if
(
kernelFunction_
->
in_support
(
dx
))
bw
=
max
(
bw
,
abs
(
j
-
map_global_to_unique
(
i
)));
}
}
}
return
bw
;
}
//-------------------------------------------------------------------
// evaluate kernels
//-------------------------------------------------------------------
void
FE_Engine
::
evaluate_kernel_functions
(
const
MATRIX
&
ptCoords
,
SPAR_MAT
&
N
)
const
{
#ifdef EXTENDED_ERROR_CHECKING
print_msg_once
(
communicator_
,
"kernel matrix bandwidth "
+
to_string
(
kernel_matrix_bandwidth
(
ptCoords
)));
#endif
if
(
!
kernelFunction_
)
throw
ATC_Error
(
"no kernel function specified"
);
int
npts
=
ptCoords
.
nRows
();
if
(
npts
>
0
)
{
N
.
reset
(
npts
,
nNodesUnique_
);
// transpose of N_Ia
DENS_VEC
xI
(
nSD_
),
x
(
nSD_
),
dx
(
nSD_
);
for
(
int
i
=
0
;
i
<
nNodes_
;
i
++
)
{
xI
=
feMesh_
->
global_coordinates
(
i
);
for
(
int
j
=
0
;
j
<
npts
;
j
++
)
{
for
(
int
k
=
0
;
k
<
nSD_
;
++
k
)
{
x
(
k
)
=
ptCoords
(
j
,
k
);
}
dx
=
x
-
xI
;
double
val
=
kernelFunction_
->
value
(
dx
);
val
*=
kernelFunction_
->
dimensionality_factor
();
if
(
val
>
0
)
N
.
add
(
j
,
map_global_to_unique
(
i
),
val
);
}
}
N
.
compress
();
}
}
//-----------------------------------------------------------------
// create a non-uniform rectangular mesh on a rectangular region
//-----------------------------------------------------------------
void
FE_Engine
::
create_mesh
(
Array
<
double
>
&
dx
,
Array
<
double
>
&
dy
,
Array
<
double
>
&
dz
,
const
char
*
regionName
,
Array
<
bool
>
periodicity
)
{
double
xmin
,
xmax
,
ymin
,
ymax
,
zmin
,
zmax
;
double
xscale
,
yscale
,
zscale
;
// check to see if region exists and is it a box, and if so get the bounds
bool
isBox
;
isBox
=
ATC
::
LammpsInterface
::
instance
()
->
region_bounds
(
regionName
,
xmin
,
xmax
,
ymin
,
ymax
,
zmin
,
zmax
,
xscale
,
yscale
,
zscale
);
if
(
!
isBox
)
throw
ATC_Error
(
"Region for FE mesh is not a box"
);
if
(
dx
(
0
)
==
0.0
)
dx
=
(
xmax
-
xmin
)
/
dx
.
size
();
if
(
dy
(
0
)
==
0.0
)
dy
=
(
ymax
-
ymin
)
/
dy
.
size
();
if
(
dz
(
0
)
==
0.0
)
dz
=
(
zmax
-
zmin
)
/
dz
.
size
();
feMesh_
=
new
FE_Rectangular3DMesh
(
dx
,
dy
,
dz
,
xmin
,
xmax
,
ymin
,
ymax
,
zmin
,
zmax
,
periodicity
,
xscale
,
yscale
,
zscale
);
stringstream
ss
;
ss
<<
"created structured mesh with "
<<
feMesh_
->
num_nodes
()
<<
" nodes, "
<<
feMesh_
->
num_nodes_unique
()
<<
" unique nodes, and "
<<
feMesh_
->
num_elements
()
<<
" elements"
;
print_msg_once
(
communicator_
,
ss
.
str
());
#ifdef ATC_VERBOSE
print_partitions
(
xmin
,
xmax
,
dx
);
print_partitions
(
ymin
,
ymax
,
dy
);
print_partitions
(
zmin
,
zmax
,
dz
);
#endif
}
//-----------------------------------------------------------------
void
FE_Engine
::
print_partitions
(
double
xmin
,
double
xmax
,
Array
<
double
>
&
dx
)
const
{
stringstream
msg
;
msg
.
precision
(
3
);
msg
<<
std
::
fixed
;
msg
<<
"
\n
index weight fraction location size[A] size[uc]:
\n
"
;
double
sum
=
0
;
for
(
int
i
=
0
;
i
<
dx
.
size
();
++
i
)
{
sum
+=
dx
(
i
);
}
double
xn
=
1.0
/
sum
;
double
xs
=
xn
*
(
xmax
-
xmin
);
double
xl
=
xs
/
(
ATC
::
LammpsInterface
::
instance
()
->
xlattice
());
double
sumn
=
0
,
sums
=
0
,
suml
=
0
;
double
x
=
xmin
;
for
(
int
i
=
0
;
i
<
dx
.
size
();
++
i
)
{
double
dxn
=
dx
(
i
)
*
xn
;
sumn
+=
dxn
;
double
dxs
=
dx
(
i
)
*
xs
;
sums
+=
dxs
;
double
dxl
=
dx
(
i
)
*
xl
;
suml
+=
dxl
;
msg
<<
std
::
setw
(
5
)
<<
i
+
1
<<
std
::
setw
(
7
)
<<
dx
(
i
)
<<
std
::
setw
(
9
)
<<
dxn
<<
std
::
setw
(
9
)
<<
x
<<
std
::
setw
(
8
)
<<
dxs
<<
std
::
setw
(
9
)
<<
dxl
<<
"
\n
"
;
x
+=
dxs
;
}
msg
<<
"sum "
<<
setw
(
7
)
<<
sum
<<
setw
(
9
)
<<
sumn
<<
setw
(
9
)
<<
x
<<
setw
(
8
)
<<
sums
<<
setw
(
9
)
<<
suml
<<
"
\n
"
;
print_msg_once
(
communicator_
,
msg
.
str
());
}
//-----------------------------------------------------------------
// create a uniform rectangular mesh on a rectangular region
//-----------------------------------------------------------------
void
FE_Engine
::
create_mesh
(
int
nx
,
int
ny
,
int
nz
,
const
char
*
regionName
,
Array
<
bool
>
periodicity
)
{
double
xmin
,
xmax
,
ymin
,
ymax
,
zmin
,
zmax
;
double
xscale
,
yscale
,
zscale
;
// check to see if region exists and is it a box, and if so get the bounds
bool
isBox
;
isBox
=
ATC
::
LammpsInterface
::
instance
()
->
region_bounds
(
regionName
,
xmin
,
xmax
,
ymin
,
ymax
,
zmin
,
zmax
,
xscale
,
yscale
,
zscale
);
if
(
!
isBox
)
throw
ATC_Error
(
"Region for FE mesh is not a box"
);
feMesh_
=
new
FE_Uniform3DMesh
(
nx
,
ny
,
nz
,
xmin
,
xmax
,
ymin
,
ymax
,
zmin
,
zmax
,
periodicity
,
xscale
,
yscale
,
zscale
);
stringstream
ss
;
ss
<<
"created uniform mesh with "
<<
feMesh_
->
num_nodes
()
<<
" nodes, "
<<
feMesh_
->
num_nodes_unique
()
<<
" unique nodes, and "
<<
feMesh_
->
num_elements
()
<<
" elements"
;
print_msg_once
(
communicator_
,
ss
.
str
());
}
//-----------------------------------------------------------------
// read a mesh from an input file using MeshReader object
//-----------------------------------------------------------------
void
FE_Engine
::
read_mesh
(
string
meshFile
,
Array
<
bool
>
&
periodicity
)
{
if
(
feMesh_
)
throw
ATC_Error
(
"FE_Engine already has a mesh"
);
feMesh_
=
MeshReader
(
meshFile
,
periodicity
).
create_mesh
();
feMesh_
->
initialize
();
}
};
Event Timeline
Log In to Comment