14 #ifndef SFCSORTBENCH_OPERATORS_H 15 #define SFCSORTBENCH_OPERATORS_H 21 #include "workspace.h" 22 #include "mathUtils.h" 23 #include "dendroProfileParams.h" 25 #define GRIDX_TO_X(xg) (((Rx/RgX)*(xg-bssn::BSSN_OCTREE_MIN[0]))+bssn::BSSN_COMPD_MIN[0]) 26 #define GRIDY_TO_Y(yg) (((Ry/RgY)*(yg-bssn::BSSN_OCTREE_MIN[1]))+bssn::BSSN_COMPD_MIN[1]) 27 #define GRIDZ_TO_Z(zg) (((Rz/RgZ)*(zg-bssn::BSSN_OCTREE_MIN[2]))+bssn::BSSN_COMPD_MIN[2]) 29 #define X_TO_GRIDX(xc) (((RgX/Rx)*(xc-bssn::BSSN_COMPD_MIN[0]))+bssn::BSSN_OCTREE_MIN[0]) 30 #define Y_TO_GRIDY(yc) (((RgY/Ry)*(yc-bssn::BSSN_COMPD_MIN[1]))+bssn::BSSN_OCTREE_MIN[1]) 31 #define Z_TO_GRIDZ(zc) (((RgZ/Rz)*(zc-bssn::BSSN_COMPD_MIN[2]))+bssn::BSSN_OCTREE_MIN[2]) 59 dendro::timer::sfcmatvec::t_elemMvec.start();
62 const double * Q1d=refEl->getQ1d();
63 const double * QT1d=refEl->getQT1d();
64 const double * Dg=refEl->getDg1d();
65 const double * W1d=refEl->getWgq();
67 const unsigned int eleOrder=refEl->getOrder();
68 const unsigned int nPe=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
69 const unsigned int nrp=eleOrder+1;
71 const double szX=fem::domain::gridX_to_X(elem.maxX())-fem::domain::gridX_to_X(elem.
minX());
72 const double szY=fem::domain::gridY_to_Y(elem.maxY())-fem::domain::gridY_to_Y(elem.minY());
73 const double szZ=fem::domain::gridZ_to_Z(elem.maxZ())-fem::domain::gridZ_to_Z(elem.minZ());
76 const double refElSz=refEl->getElementSz();
79 DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,Q1d,f_rhs,imV1);
80 DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,Q1d,imV1,imV2);
81 DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,Q1d,imV2,out);
85 const double Jx = 1.0/(refElSz/(double (szX)));
86 const double Jy = 1.0/(refElSz/(double (szY)));
87 const double Jz = 1.0/(refElSz/(double (szZ)));
91 for(
unsigned int k=0;k<(eleOrder+1);k++)
92 for(
unsigned int j=0;j<(eleOrder+1);j++)
93 for(
unsigned int i=0;i<(eleOrder+1);i++)
94 out[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]*=(Jx*Jy*Jz*W1d[i]*W1d[j]*W1d[k]);
98 DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,QT1d,out,imV1);
99 DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,QT1d,imV1,imV2);
100 DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,QT1d,imV2,out);
102 #ifdef FEM_ACCUMILATE_ONES_TEST 103 for(
unsigned int k=0;k<(eleOrder+1);k++)
104 for(
unsigned int j=0;j<(eleOrder+1);j++)
105 for(
unsigned int i=0;i<(eleOrder+1);i++)
106 out[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]=1.0;
119 #ifdef MATVEC_PROFILE 120 dendro::timer::sfcmatvec::t_elemMvec.stop();
129 template <
typename T>
133 #ifdef MATVEC_PROFILE 134 dendro::timer::sfcmatvec::t_elemMvec.start();
137 const double * Q1d=refEl->getQ1d();
138 const double * QT1d=refEl->getQT1d();
139 const double * Dg=refEl->getDg1d();
140 const double * DgT=refEl->getDgT1d();
141 const double * W1d=refEl->getWgq();
143 const unsigned int eleOrder=refEl->getOrder();
144 const unsigned int nPe=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
145 const unsigned int nrp=eleOrder+1;
147 const unsigned int sz=1u<<(m_uiMaxDepth-elem.
getLevel());
148 const double refElSz=refEl->getElementSz();
150 DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,Dg,in,imV1);
151 DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,Q1d,imV1,imV2);
152 DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,Q1d,imV2,Qx);
155 DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,Q1d,in,imV1);
156 DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,Dg,imV1,imV2);
157 DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,Q1d,imV2,Qy);
160 DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,Q1d,in,imV1);
161 DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,Q1d,imV1,imV2);
162 DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,Dg,imV2,Qz);
165 const double szX=fem::domain::gridX_to_X(elem.maxX())-fem::domain::gridX_to_X(elem.
minX());
166 const double szY=fem::domain::gridY_to_Y(elem.maxY())-fem::domain::gridY_to_Y(elem.minY());
167 const double szZ=fem::domain::gridZ_to_Z(elem.maxZ())-fem::domain::gridZ_to_Z(elem.minZ());
169 const double Jx = 1.0/(refElSz/(double (szX)));
170 const double Jy = 1.0/(refElSz/(double (szY)));
171 const double Jz = 1.0/(refElSz/(double (szZ)));
175 for(
unsigned int k=0;k<(eleOrder+1);k++)
176 for(
unsigned int j=0;j<(eleOrder+1);j++)
177 for(
unsigned int i=0;i<(eleOrder+1);i++)
179 Qx[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]*=( ((Jy*Jz)/Jx)*W1d[i]*W1d[j]*W1d[k]);
180 Qy[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]*=( ((Jx*Jz)/Jy)*W1d[i]*W1d[j]*W1d[k]);
181 Qz[k*(eleOrder+1)*(eleOrder+1)+j*(eleOrder+1)+i]*=( ((Jx*Jy)/Jz)*W1d[i]*W1d[j]*W1d[k]);
186 DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,DgT,Qx,imV1);
187 DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,QT1d,imV1,imV2);
188 DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,QT1d,imV2,Qx);
190 DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,QT1d,Qy,imV1);
191 DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,DgT,imV1,imV2);
192 DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,QT1d,imV2,Qy);
194 DENDRO_TENSOR_IIAX_APPLY_ELEM(nrp,QT1d,Qz,imV1);
195 DENDRO_TENSOR_IAIX_APPLY_ELEM(nrp,QT1d,imV1,imV2);
196 DENDRO_TENSOR_AIIX_APPLY_ELEM(nrp,DgT,imV2,Qz);
198 for(
unsigned int i=0;i<nPe;i++)
199 out[i]=Qx[i]+Qy[i]+Qz[i];
202 #ifdef MATVEC_PROFILE 203 dendro::timer::sfcmatvec::t_elemMvec.stop();
211 template <
typename T>
212 void applyDirichletBoundary(
const ot::Mesh* mesh,T* in,T bdyValue)
219 const unsigned int * e2n=&(*(mesh->
getE2NMapping().begin()));
222 const unsigned int d_min=0;
223 const unsigned int d_max=(1u<<m_uiMaxDepth);
227 const unsigned int nx=(eleOrder+1);
228 const unsigned int ny=(eleOrder+1)*(eleOrder+1);
229 const unsigned int nz=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
230 const double bdy_value=bdyValue;
235 x=pNodes[ele].
minX();
236 y=pNodes[ele].minY();
237 z=pNodes[ele].minZ();
241 for(
unsigned int k=0;k<(eleOrder+1);k++)
242 for(
unsigned int j=0;j<(eleOrder+1);j++)
245 in[e2n[ele*nPe+k*ny+j*nx+0]]=bdy_value;
251 for(
unsigned int k=0;k<(eleOrder+1);k++)
252 for(
unsigned int i=0;i<(eleOrder+1);i++)
255 in[e2n[ele*nPe+k*ny+0*nx+i]]=bdy_value;
262 for(
unsigned int j=0;j<(eleOrder+1);j++)
263 for(
unsigned int i=0;i<(eleOrder+1);i++)
266 in[e2n[ele*nPe+0*ny+j*nx+i]]=bdy_value;
270 x=pNodes[ele].maxX();
271 y=pNodes[ele].maxY();
272 z=pNodes[ele].maxZ();
276 for(
unsigned int k=0;k<(eleOrder+1);k++)
277 for(
unsigned int j=0;j<(eleOrder+1);j++)
280 in[e2n[ele*nPe+k*ny+j*nx+eleOrder]]=bdy_value;
286 for(
unsigned int k=0;k<(eleOrder+1);k++)
287 for(
unsigned int i=0;i<(eleOrder+1);i++)
290 in[e2n[ele*nPe+k*ny+eleOrder*nx+i]]=bdy_value;
297 for(
unsigned int j=0;j<(eleOrder+1);j++)
298 for(
unsigned int i=0;i<(eleOrder+1);i++)
301 in[e2n[ele*nPe+eleOrder*ny+j*nx+i]]=bdy_value;
320 template <
typename T>
321 void matvec(
ot::Mesh* mesh,T* in, T* out,
bool applyBoundary=
true)
334 T* elementalVec=
new T[nPe];
335 T* elementalMVec=
new T[nPe];
337 const unsigned int nx=(eleOrder+1);
338 const unsigned int ny=(eleOrder+1)*(eleOrder+1);
339 const unsigned int nz=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
343 unsigned int nodeLookup,ownerID,ii_x,jj_y,kk_z;
350 for(
unsigned int gIndex=0;gIndex<level1Ghost.size();gIndex++)
352 unsigned int ele=level1Ghost[gIndex];
354 operators::poisson::elementStiffnessMatvec(elementalVec,pNodes[ele],refEl,elementalMVec);
369 operators::poisson::elementStiffnessMatvec(elementalVec,pNodes[ele],refEl,elementalMVec);
379 delete [] elementalVec;
380 delete [] elementalMVec;
395 template <
typename T>
396 void computeRHS(
ot::Mesh* mesh,T* in, T* out)
409 T* elementalVec=
new T[nPe];
410 T* elementalMVec=
new T[nPe];
412 const unsigned int nx=(eleOrder+1);
413 const unsigned int ny=(eleOrder+1)*(eleOrder+1);
414 const unsigned int nz=(eleOrder+1)*(eleOrder+1)*(eleOrder+1);
422 const unsigned int * e2n=&(*(mesh->
getE2NMapping().begin()));
424 unsigned int nodeLookup,ownerID,ii_x,jj_y,kk_z;
433 for(
unsigned int gIndex=0;gIndex<level1Ghost.size();gIndex++)
435 unsigned int ele=level1Ghost[gIndex];
437 #ifdef MATVEC_PROFILE 438 dendro::timer::sfcmatvec::t_p2cInterp.start();
441 #ifdef MATVEC_PROFILE 442 dendro::timer::sfcmatvec::t_p2cInterp.stop();
445 operators::poisson::elementMassMatvec(elementalVec,pNodes[ele],refEl,elementalMVec);
448 #ifdef MATVEC_PROFILE 449 dendro::timer::sfcmatvec::t_c2pInterp.start();
454 #ifdef MATVEC_PROFILE 455 dendro::timer::sfcmatvec::t_c2pInterp.stop();
464 #ifdef MATVEC_PROFILE 465 dendro::timer::sfcmatvec::t_p2cInterp.start();
469 #ifdef MATVEC_PROFILE 470 dendro::timer::sfcmatvec::t_p2cInterp.stop();
472 operators::poisson::elementMassMatvec(elementalVec,pNodes[ele],refEl,elementalMVec);
475 #ifdef MATVEC_PROFILE 476 dendro::timer::sfcmatvec::t_c2pInterp.start();
480 #ifdef MATVEC_PROFILE 481 dendro::timer::sfcmatvec::t_c2pInterp.stop();
489 delete [] elementalVec;
490 delete [] elementalMVec;
508 #endif //SFCSORTBENCH_OPERATORS_H const std::vector< unsigned int > & getLevel1GhostElementIndices() const
returns the Level 1 ghost element indices.
Definition: mesh.h:1061
void computeElementalContribution(const T *in, T *out, unsigned int elementID) const
: Computes the contribution of elemental nodal values to the parent elements if it is hanging...
unsigned int minX() const
returns true min corner(anchor) and the max corner coordinates of a octant.
Definition: TreeNode.h:164
A class to manage octants.
Definition: TreeNode.h:35
void performGhostExchange(std::vector< T > &vec)
Perform the ghost exchange for the vector vec.
const std::vector< unsigned int > & getE2NMapping_DG() const
Definition: mesh.h:1082
unsigned int getNodeLocalBegin() const
return the location of local node begin
Definition: mesh.h:1043
bool isNodeLocal(unsigned int eleID, unsigned int ix, unsigned int jy, unsigned int kz) const
: Returns true if the specified node (e,i,j,k) is local.
Definition: mesh.h:1358
workspace variables allocations
Definition: matvec.h:44
unsigned int getElementLocalBegin() const
return the begin location of element local
Definition: mesh.h:1030
void getElementNodalValues(const T *vec, T *nodalValues, unsigned int elementID) const
: Returns the nodal values of a given element for a given variable vector.
const std::vector< ot::TreeNode > & getAllElements() const
returns the pointer to All elements array.
Definition: mesh.h:1058
unsigned int getNumNodesPerElement() const
return the number of nodes per element.
Definition: mesh.h:1090
const RefElement * getReferenceElement() const
returns const pointer to reference element
Definition: mesh.h:1114
unsigned int getLevel() const
return the level of the of the octant
Definition: TreeNode.h:387
bool isNodeHanging(unsigned int eleID, unsigned int ix, unsigned int jy, unsigned int kz) const
: Returns true if the specified node (e,i,j,k) is hanging.
Definition: mesh.cpp:8943
const std::vector< unsigned int > & getE2NMapping() const
Definition: mesh.h:1073
bool isActive() const
returns true if mesh is active
Definition: mesh.h:1173
unsigned int getNodeLocalEnd() const
return the location of local node end
Definition: mesh.h:1045
unsigned int getElementOrder() const
returns the order of an element
Definition: mesh.h:1093
unsigned int getElementLocalEnd() const
return the end location of element local
Definition: mesh.h:1032