27 #ifndef _SF_CONTAINER_H
28 #define _SF_CONTAINER_H
74 return {arr[0], arr[1], arr[2]};
80 arr[0] = p.
x, arr[1] = p.
y, arr[2] = p.
z;
92 return a.
x*b.
x + a.
y*b.
y + a.
z*b.
z;
95 inline void outer_prod(
const Point & a,
const Point & b,
const double s,
double* buff,
const bool add =
false)
98 memset(buff, 0, 9*
sizeof(
double));
100 buff[0] += a.
x*b.
x*s, buff[1] += a.
x*b.
y*s, buff[2] += a.
x*b.
z*s;
101 buff[3] += a.
y*b.
x*s, buff[4] += a.
y*b.
y*s, buff[5] += a.
y*b.
z*s;
102 buff[6] += a.
z*b.
x*s, buff[7] += a.
z*b.
y*s, buff[8] += a.
z*b.
z*s;
118 double norm =
mag(vect);
119 return {vect.
x / norm, vect.
y / norm, vect.
z / norm};
124 return {a.
x - b.
x, a.
y - b.
y, a.
z - b.
z};
129 return {a.
x + b.
x, a.
y + b.
y, a.
z + b.
z};
134 return {a.
x * s, a.
y * s, a.
z * s};
159 if ( !strcmp( eletype,
"Tt" ) ) {
161 }
else if ( !strcmp( eletype,
"Hx" ) ) {
163 }
else if ( !strcmp( eletype,
"Oc" ) ) {
165 }
else if ( !strcmp( eletype,
"Py" ) ) {
167 }
else if ( !strcmp( eletype,
"Pr" ) ) {
169 }
else if ( !strcmp( eletype,
"Qd" ) ) {
171 }
else if ( !strcmp( eletype,
"Tr" ) ) {
179 template<
class S,
class POINT>
inline
183 for(
size_t i=0; i<pts.
size(); i++)
184 pts[i] = {arr[i*3+0], arr[i*3+1], arr[i*3+2]};
197 template<
class T>
class overlapping_layout;
198 template<
class T>
class non_overlapping_layout;
243 for(
size_t i=0; i<a.
size(); i++)
245 T aidx = a[i], bidx = b[i];
282 size_t lsize = idx.
size();
296 T val = it != _fwd.
end() ? it->second : T(-1);
298 while( (c < lsize) && (cc == idx[c]) )
318 size_t lsize = idx.
size();
332 T val = it != _bwd.
end() ? *it : T(-1);
334 while( (c < lsize) && (cc == idx[c]) )
351 return _fwd.
count(idx) == 1;
355 return _bwd.
count(idx) == 1;
382 template<
class T,
class S>
414 std::map< SF_nbr, vector<T> >
nbr;
434 if(
nbr.count(nbr_type) == 0)
437 return nbr[nbr_type];
454 typename std::map< SF_nbr, vector<T> >::iterator it =
nbr.find(nbr_type);
455 assert(it !=
nbr.end());
472 typename std::map< SF_nbr, vector<T> >::const_iterator it =
nbr.find(nbr_type);
473 assert(it !=
nbr.end());
490 for(
const T & c : this->con) g2l[c] = 0;
497 for(
auto it = g2l.
begin(); it != g2l.
end(); ++it, widx++) {
498 nod[widx] = it->first;
503 this->l_numpts = nod.
size();
517 for(T & c :
con) c = nod[c];
547 con.resize(0);
con.reallocate();
548 dsp.resize(0);
dsp.reallocate();
549 tag.resize(0);
tag.reallocate();
570 template<
class T,
class S>
595 template<
class T,
class S>
602 for(
size_t nidx=0; nidx < n2n_cnt.
size(); nidx++)
603 if(nmax < n2n_cnt[nidx]) nmax = n2n_cnt[nidx];
605 MPI_Allreduce(MPI_IN_PLACE, &nmax, 1, MPI_INT, MPI_MAX, mesh.
comm);
610 template<
class T,
class S>
613 const vector<T> & alg_nod = mesh.
pl.algebraic_nodes();
616 for(
const T & n : alg_nod) alg_mask[n] =
true;
636 scnt.assign(size, T()),
rcnt.assign(size, T()),
sdsp.assign(size+1, T()),
rdsp.assign(size+1, T());
643 for(
size_t i=0; i<
scnt.size(); i++)
646 for(
size_t i=0; i<
sdsp.size(); i++)
670 MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
675 MPI_Alltoall(
scnt.data(),
sizeof(T), MPI_BYTE,
rcnt.data(),
sizeof(T), MPI_BYTE, comm);
692 MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
697 for(
size_t i=0; i<dest.
size(); i++)
scnt[dest[i]] += cnt[i];
699 MPI_Alltoall(
scnt.data(),
sizeof(T), MPI_BYTE,
rcnt.data(),
sizeof(T), MPI_BYTE, comm);
716 size_t size =
rcnt.size();
719 for(
size_t i=0, widx=0; i<size; i++)
720 for(T j=0; j<
rcnt[i]; j++, widx++)
731 template<
class T,
class S>
745 template<
class T,
class S,
class V>
764 return lhs.
v1 < rhs.
v1;
766 return lhs.
v2 < rhs.
v2;
768 template<
class T,
class S>
772 return lhs.
v1 < rhs.
v1;
774 return lhs.
v2 < rhs.
v2;
781 return lhs.
v1 < rhs.
v1;
782 else if(lhs.
v2 != rhs.
v2)
783 return lhs.
v2 < rhs.
v2;
785 return lhs.
v3 < rhs.
v3;
792 return lhs.
v1 < rhs.
v1;
793 else if(lhs.
v2 != rhs.
v2)
794 return lhs.
v2 < rhs.
v2;
795 else if(lhs.
v3 != rhs.
v3)
796 return lhs.
v3 < rhs.
v3;
798 return lhs.
v4 < rhs.
v4;
812 return a.
v1 == b.
v1 && a.
v2 == b.
v2;
881 void sort_triple(
const T in1,
const T in2,
const T in3, T & out1, T & out2, T & out3)
883 bool t12 = in1 < in2;
884 bool t13 = in1 < in3;
885 bool t23 = in2 < in3;
891 out1=in1; out2=in2; out3=in3;
897 out1=in1; out2=in3; out3=in2;
901 out1=in3; out2=in1; out3=in2;
911 out1=in2; out2=in1; out3=in3;
915 out1=in2; out2=in3; out3=in1;
920 out1=in3; out2=in2; out3=in1;
void sort_triple(const T in1, const T in2, const T in3, T &out1, T &out2, T &out3)
sort the "in" triple into the "out" triple
#define SF_COMM
the default SlimFem MPI communicator
Sequential linear algebra kernels.
The vector class and related algorithms.
The class holds the communication graph for a MPI_Exchange() call.
vector< T > rcnt
Number of elements received from each rank.
vector< T > scnt
Number of elements sent to each rank.
void resize(size_t size)
Resize all vectors to size.
void configure(const vector< V > &dest, MPI_Comm comm)
Set up the communication graph.
void configure(const vector< V > &dest, const vector< T > &cnt, MPI_Comm comm)
Set up the communication graph.
vector< T > sdsp
Displacements w.r.t. scnt.
vector< T > rdsp
Displacements w.r.t. rcnt.
void source_ranks(vector< V > &source)
For every received data element, get the rank indices it was receive from.
void transpose()
transpose comm graph (receive becomes send, and vice versa)
void scale(V fac)
scale comm graph layout data
Index mapping class. This is a bijective mapping.
void assign(const vector< T > &a, const vector< T > &b)
Set up the index mapping between a and b.
const hashmap::unordered_map< T, T > & get_fwd_map() const
void backward_map(vector< T > &idx) const
Map a whole array of indices in b to indices in a.
hashmap::unordered_map< T, T > & get_bwd_map()
index_mapping(const vector< T > &a, const vector< T > &b)
Constructor that uses assign() to set up the index mapping.
size_t size() const
number of entries in bijective map
void forward_map(vector< T > &idx) const
Map a whole array of indices in a to indices in b.
hashmap::unordered_map< T, T > & get_fwd_map()
const hashmap::unordered_map< T, T > & get_bwd_map() const
T forward_map(T idx) const
Map one index from a to b.
bool in_a(const T idx)
return whether idx is in set A
bool in_b(const T idx)
return whether idx is in set B
T backward_map(T idx) const
Map one index from b to a.
index_mapping()
empty constructor.
The mesh storage class. It contains both element and vertex data.
void globalize(SF_nbr nbr_type)
Localize the connectivity data w.r.t. a given numbering.
void clear_data()
Clear the mesh data from memory.
overlapping_layout< T > pl
nodal parallel layout
size_t g_numpts
global number of points
meshdata()
construct empty mesh
vector< T > dsp
connectivity starting index of each element
vector< S > she
sheet direction
vector< S > fib
fiber direction
size_t l_numelem
local number of elements
vector< elem_t > type
element type
std::map< SF_nbr, vector< T > > nbr
container for different numberings
std::string name
the mesh name
void generate_par_layout()
Set up the parallel layout.
vector< T > & register_numbering(SF_nbr nbr_type)
Register a new numbering to the mesh and return the associated index vector.
vector< S > xyz
node cooridnates
const vector< T > & get_numbering(SF_nbr nbr_type) const
Get the vector defining a certain numbering.
size_t l_numpts
local number of points
size_t g_numelem
global number of elements
void localize(SF_nbr nbr_type)
Localize the connectivity data w.r.t. a given numbering.
MPI_Comm comm
the parallel mesh is defined on a MPI world
vector< T > & get_numbering(SF_nbr nbr_type)
Get the vector defining a certain numbering.
vector< T > tag
element tag
hashmap::unordered_set< int > extr_tag
the element tags based on which the mesh has been extracted
non_overlapping_layout< T > epl
element parallel layout
The parallel layout of non overlapping indices.
The overlapping_layout class contains the algorithms related to managing overlapping parallel index s...
void assign(const vector< T > &idx, MPI_Comm comm)
Initialization function.
size_t size() const
The current size of the vector.
void resize(size_t n)
Resize a vector.
void assign(InputIterator s, InputIterator e)
Assign a memory range.
iterator find(const K &key)
Search for key. Return iterator.
hm_int count(const K &key) const
Check if key exists.
void sort(Compare comp=Compare())
Sort data entries.
Dense matrix class and associated funcs.
Classes similar to unordered_set and unordered_map, but with better performance.
unsigned long int hm_uint
double mag(const Point &vect)
vector magnitude
void cnt_from_dsp(const vector< T > &dsp, vector< T > &cnt)
Compute counts from displacements.
Point arr_to_point(V *arr)
void dsp_from_cnt(const vector< T > &cnt, vector< T > &dsp)
Compute displacements from counts.
void transpose_connectivity(const vector< T > &a_cnt, const vector< T > &a_con, vector< T > &b_cnt, vector< T > &b_con)
Transpose CRS matrix graph A into B.
bool operator<(const struct tuple< T > &lhs, const struct tuple< T > &rhs)
dmat< S > operator-(const dmat< S > &a, const dmat< S > &b)
void interval(vector< T > &vec, size_t start, size_t end)
Create an integer interval between start and end.
void get_alg_mask(const meshdata< T, S > &mesh, vector< bool > &alg_mask)
double inner_prod(const Point &a, const Point &b)
void binary_sort_copy(vector< T > &_V, vector< S > &_W)
T sum(const vector< T > &vec)
Compute sum of a vector's entries.
void count(const vector< T > &data, vector< S > &cnt)
Count number of occurrences of indices.
void point_to_arr(Point &p, V *arr)
Point project(const Point &a, const Point &b)
project b onto a
void outer_prod(const Point &a, const Point &b, const double s, double *buff, const bool add=false)
void multiply_connectivities(const vector< T > &a_cnt, const vector< T > &a_con, const vector< T > &b_cnt, const vector< T > &b_con, vector< T > &c_cnt, vector< T > &c_con)
int max_nodal_edgecount(const meshdata< T, S > &mesh)
Compute the maximum number of node-to-node edges for a mesh.
dmat< S > operator*(const dmat< S > &a, const dmat< S > &b)
Point normalize(const Point &vect)
void array_to_points(const vector< S > &arr, vector< POINT > &pts)
Point orthogonalize(const Point &a, const Point &b)
void global_to_local(const vector< T > &glob, vector< T > &data, bool sortedData, bool doWarn)
elem_t getElemTypeID(char *eletype)
Generate element type enum from string.
double distance(const Point &a, const Point &b)
dmat< S > operator+(const dmat< S > &a, const dmat< S > &b)
Point cross(const Point &a, const Point &b)
cross product
void nodal_connectivity_graph(const meshdata< T, S > &mesh, vector< T > &n2n_cnt, vector< T > &n2n_con)
Compute the node-to-node connectivity.
SF_nbr
Enumeration encoding the different supported numberings.
@ NBR_PETSC
PETSc numbering of nodes.
@ NBR_ELEM_REF
The element numbering of the reference mesh (the one stored on HD).
@ NBR_REF
The nodal numbering of the reference mesh (the one stored on HD).
@ NBR_SUBMESH
Submesh nodal numbering: The globally ascending sorted reference indices are reindexed.
@ NBR_ELEM_SUBMESH
Submesh element numbering: The globally ascending sorted reference indices are reindexed.
const hm_uint mkhash_init
hm_uint mkhash(hm_uint a, hm_uint b)
static hm_uint hash(SF::quadruple< T > a)
static bool cmp(SF::quadruple< T > a, SF::quadruple< T > b)
static bool cmp(SF::triple< T > a, SF::triple< T > b)
static hm_uint hash(SF::triple< T > a)
static hm_uint hash(SF::tuple< T > a)
static bool cmp(SF::tuple< T > a, SF::tuple< T > b)