50 namespace MultiRegions
77 ExpList1D::ExpList1D():
88 ExpList(In,DeclareCoeffPhysArrays)
97 const std::vector<unsigned int> &eIDs,
98 const bool DeclareCoeffPhysArrays,
100 ExpList(In, eIDs, DeclareCoeffPhysArrays)
145 = graph1D->GetExpansions();
149 SpatialDomains::ExpansionMap::const_iterator expIt;
150 for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
152 if ((SegmentGeom = boost::dynamic_pointer_cast<
154 expIt->second->m_geomShPtr)))
158 seg->SetElmtId(
id++);
159 (*m_exp).push_back(seg);
163 ASSERTL0(
false,
"dynamic cast to a SegGeom failed");
205 const bool DeclareCoeffPhysArrays,
217 = graph1D->GetExpansions();
220 SpatialDomains::ExpansionMap::const_iterator expIt;
221 for (expIt = expansions.begin(); expIt != expansions.end(); ++expIt)
226 if ((SegmentGeom = boost::dynamic_pointer_cast<
228 expIt->second->m_geomShPtr)))
234 seg->SetElmtId(
id++);
237 (*m_exp).push_back(seg);
241 ASSERTL0(
false,
"dynamic cast to a SegGeom failed");
253 if(DeclareCoeffPhysArrays)
290 const bool DeclareCoeffPhysArrays,
291 const std::string var,
292 bool SetToOneSpaceDimension,
300 SpatialDomains::CompositeMap::const_iterator compIt;
304 = graph1D->GetExpansions(var);
307 for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
309 comp = compIt->second;
312 for(j = 0; j < compIt->second->size(); ++j)
314 SpatialDomains::ExpansionMap::const_iterator expIt;
316 if((SegmentGeom = boost::dynamic_pointer_cast<
318 (*compIt->second)[j])))
321 if((expIt = expansions.find(SegmentGeom->GetGlobalID())) != expansions.end())
325 if(SetToOneSpaceDimension)
328 SegmentGeom->GenerateOneSpaceDimGeom();
341 ASSERTL0(
false,
"Failed to find basis key");
346 ASSERTL0(
false,
"Failed to dynamic cast geometry to SegGeom");
350 seg->SetElmtId(
id++);
353 (*m_exp).push_back(seg);
365 if(DeclareCoeffPhysArrays)
393 const bool DeclareCoeffPhysArrays,
394 const std::string variable,
404 SpatialDomains::CompositeMap::const_iterator compIt;
409 for(compIt = domain.begin(); compIt != domain.end(); ++compIt)
411 comp = compIt->second;
413 for(j = 0; j < compIt->second->size(); ++j)
415 if((SegmentGeom = boost::dynamic_pointer_cast<
417 (*compIt->second)[j])))
427 seg->SetElmtId(
id++);
428 (*m_exp).push_back(seg);
432 ASSERTL0(
false,
"dynamic cast to a SegGeom failed");
447 if(DeclareCoeffPhysArrays)
479 const bool DeclareCoeffPhysArrays,
480 const std::string variable,
484 int i, j, id, elmtid = 0;
496 map<int,int> EdgeDone;
497 map<int,int> NormalSet;
503 for(i = 0; i < bndCond.num_elements(); ++i)
505 if(bndCond[i]->GetBoundaryConditionType()
508 for(j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
511 ->GetExp(j)->GetBasis(0)->GetBasisKey();
512 exp1D = bndConstraint[i]->GetExp(j)->
513 as<LocalRegions::Expansion1D>();
514 segGeom = exp1D->GetGeom1D();
518 edgesDone.insert(segGeom->GetEid());
520 seg->SetElmtId(elmtid++);
521 (*m_exp).push_back(seg);
531 for(i = 0; i < locexp.size(); ++i)
535 for(j = 0; j < locexp[i]->GetNedges(); ++j)
537 segGeom = exp2D->
GetGeom2D()->GetEdge(j);
538 id = segGeom->GetEid();
540 if (edgesDone.count(
id) != 0)
545 it = edgeOrders.find(
id);
547 if (it == edgeOrders.end())
549 edgeOrders.insert(std::make_pair(
id, std::make_pair(
550 segGeom, locexp[i]->DetEdgeBasisKey(j))));
555 = locexp[i]->DetEdgeBasisKey(j);
564 if (np2 >= np1 && nm2 >= nm1)
568 else if (np2 < np1 && nm2 < nm1)
570 it->second.second = edge;
575 "inappropriate number of points/modes (max "
576 "num of points is not set with max order)");
583 pSession->GetComm()->GetRowComm();
584 int nproc = vComm->GetSize();
585 int edgepr = vComm->GetRank();
592 for(i = 0; i < locexp.size(); ++i)
594 eCnt += locexp[i]->GetNedges();
600 edgesCnt[edgepr] = eCnt;
606 for (i = 1; i < nproc; ++i)
608 eTotOffsets[i] = eTotOffsets[i-1] + edgesCnt[i-1];
616 int cntr = eTotOffsets[edgepr];
618 for(i = 0; i < locexp.size(); ++i)
624 for(j = 0; j < nedges; ++j, ++cntr)
627 locexp[i]->DetEdgeBasisKey(j);
628 EdgesTotID [cntr] = exp2D->GetGeom2D()->GetEid(j);
638 for (i = 0; i < totEdgeCnt; ++i)
640 it = edgeOrders.find(EdgesTotID[i]);
642 if (it == edgeOrders.end())
660 if (np2 >= np1 && nm2 >= nm1)
664 else if (np2 < np1 && nm2 < nm1)
666 it->second.second = edge;
671 "inappropriate number of points/modes (max "
672 "num of points is not set with max order)");
677 for (it = edgeOrders.begin(); it != edgeOrders.end(); ++it)
681 seg->SetElmtId(elmtid++);
682 (*m_exp).push_back(seg);
694 if(DeclareCoeffPhysArrays)
733 int quad_npoints = elmExp->GetTotPoints();
735 elmExp->GetPointsType(0));
742 int kernel_width = kernel->GetKernelWidth();
749 for(i = 0; i < inarray.num_elements(); i++)
752 kernel->MoveKernelCenter(inarray[i],local_kernel_breaks);
756 kernel->FindMeshUnderKernel(local_kernel_breaks,h,mesh_breaks);
759 int total_nbreaks = local_kernel_breaks.num_elements() +
760 mesh_breaks.num_elements();
763 kernel->Sort(local_kernel_breaks,mesh_breaks,total_breaks);
768 for(j = 0; j < total_breaks.num_elements()-1; j++)
774 for(r = 0; r < quad_points.num_elements(); r++)
776 mapped_quad_points[r]
777 = (quad_points[r] + 1.0) * 0.5 * (b - a) + a;
786 elmExp->GetBasisNumModes(0),u_value);
790 kernel->EvaluateKernel(mapped_quad_points,h,k_value);
793 for(r = 0; r < quad_npoints; r++)
795 integral_value += (b - a) * 0.5 * k_value[r]
796 * u_value[r] * quad_weights[r];
799 outarray[i] = integral_value/h;
830 for(i = 0; i < outarray.num_elements(); i++)
836 int x_size = inarray2.num_elements();
841 for(i = 0; i < x_size; i++ )
843 x_elm[i] = (int)floor(inarray2[i]/h);
847 for(i = 0; i < x_size; i++)
853 while(x_elm[i] >= num_elm)
855 x_elm[i] -= num_elm ;
860 for(i = 0; i < x_size; i++)
862 x_values_cp[i] = (inarray2[i]/h - floor(inarray2[i]/h))*2 - 1.0;
869 for(i = 0; i < nmodes; i++)
872 jacobi_poly.get()+i*x_size,NULL,i,0.0,0.0);
876 for(r = 0; r < nmodes; r++)
878 for(j = 0; j < x_size; j++)
880 int index = ((x_elm[j])*nmodes)+r;
881 outarray[j] += inarray1[index]*jacobi_poly[r][j];
925 for (i = 0; i <
m_exp->size(); ++i)
927 for (j = 0; j < (*m_exp)[i]->GetNverts(); ++j)
929 (*m_exp)[i]->ComputeVertexNormal(j);
950 int i,j,k,e_npoints,offset;
957 ASSERTL1(Vec.num_elements() >= coordim,
958 "Input vector does not have sufficient dimensions to "
962 for(i = 0; i <
m_exp->size(); ++i)
965 e_npoints = (*m_exp)[i]->GetNumPoints(0);
966 normals = (*m_exp)[i]->GetPhysNormals();
972 for(j = 0; j < e_npoints; ++j)
976 for(k = 0; k < coordim; ++k)
978 Vn += Vec[k][offset+j]*normals[k*e_npoints + j];
984 Upwind[offset + j] = Fwd[offset + j];
988 Upwind[offset + j] = Bwd[offset + j];
1013 int i,j,e_npoints,offset;
1017 for(i = 0; i <
m_exp->size(); ++i)
1020 e_npoints = (*m_exp)[i]->GetNumPoints(0);
1024 for(j = 0; j < e_npoints; ++j)
1027 if(Vn[offset + j] > 0.0)
1029 Upwind[offset + j] = Fwd[offset + j];
1033 Upwind[offset + j] = Bwd[offset + j];
1050 int i,j,k,e_npoints,offset;
1058 ASSERTL1(normals.num_elements() >= coordim,
1059 "Output vector does not have sufficient dimensions to "
1062 for (i = 0; i <
m_exp->size(); ++i)
1069 int edgeNumber = loc_exp->GetLeftAdjacentElementEdge();
1072 e_npoints = (*m_exp)[i]->GetNumPoints(0);
1074 locnormals = loc_elmt->GetEdgeNormal(edgeNumber);
1075 int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
1076 int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
1078 if (e_nmodes != loc_nmodes)
1080 if (loc_exp->GetRightAdjacentElementEdge() >= 0)
1083 loc_exp->GetRightAdjacentElementExp();
1085 int EdgeNumber = loc_exp->GetRightAdjacentElementEdge();
1088 locnormals = loc_elmt->GetEdgeNormal(EdgeNumber);
1093 for (j = 0; j < e_npoints; ++j)
1097 for (k = 0; k < coordim; ++k)
1099 normals[k][offset + j] = -locnormals[k][j];
1108 for (
int p = 0;
p < coordim; ++
p)
1112 loc_exp->GetBasis(0)->GetPointsKey();
1114 loc_elmt->GetBasis(0)->GetPointsKey();
1124 for (j = 0; j < e_npoints; ++j)
1128 for (k = 0; k < coordim; ++k)
1130 normals[k][offset + j] = normal[k][j];
1141 for (j = 0; j < e_npoints; ++j)
1145 for (k = 0; k < coordim; ++k)
1147 normals[k][offset + j] = locnormals[k][j];
1179 int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
1181 int ntotminus = (nquad0-1);
1187 (*m_exp)[expansion]->GetCoords(coords[0],coords[1],coords[2]);
1189 outfile <<
" <Piece NumberOfPoints=\""
1190 << ntot <<
"\" NumberOfCells=\""
1191 << ntotminus <<
"\">" << endl;
1192 outfile <<
" <Points>" << endl;
1193 outfile <<
" <DataArray type=\"Float64\" "
1194 <<
"NumberOfComponents=\"3\" format=\"ascii\">" << endl;
1196 for (i = 0; i < ntot; ++i)
1198 for (j = 0; j < 3; ++j)
1200 outfile << setprecision(8) << scientific
1201 << (float)coords[j][i] <<
" ";
1206 outfile <<
" </DataArray>" << endl;
1207 outfile <<
" </Points>" << endl;
1208 outfile <<
" <Cells>" << endl;
1209 outfile <<
" <DataArray type=\"Int32\" "
1210 <<
"Name=\"connectivity\" format=\"ascii\">" << endl;
1211 for (i = 0; i < nquad0-1; ++i)
1213 outfile << i <<
" " << i+1 << endl;
1216 outfile <<
" </DataArray>" << endl;
1217 outfile <<
" <DataArray type=\"Int32\" "
1218 <<
"Name=\"offsets\" format=\"ascii\">" << endl;
1219 for (i = 0; i < ntotminus; ++i)
1221 outfile << i*2+2 <<
" ";
1224 outfile <<
" </DataArray>" << endl;
1225 outfile <<
" <DataArray type=\"UInt8\" "
1226 <<
"Name=\"types\" format=\"ascii\">" << endl;
1227 for (i = 0; i < ntotminus; ++i)
1232 outfile <<
" </DataArray>" << endl;
1233 outfile <<
" </Cells>" << endl;
1234 outfile <<
" <PointData>" << endl;
const Array< OneD, const NekDouble > & GetCoeffs() const
This function returns (a reference to) the array (implemented as m_coeffs) containing all local expa...
ExpList1D()
The default constructor.
#define ASSERTL0(condition, msg)
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
NekOptimize::GlobalOptParamSharedPtr m_globalOptParam
BasisType GetBasisType() const
Return type of expansion basis.
const boost::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
PointsType GetPointsType() const
Return type of quadrature.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
int GetExpSize(void)
This function returns the number of elements in the expansion.
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
virtual void v_SetUpPhysNormals()
Set up the normals on each expansion.
std::vector< ExpansionSharedPtr > ExpansionVector
Expansion2DSharedPtr GetLeftAdjacentElementExp() const
void v_Upwind(const Array< OneD, const Array< OneD, NekDouble > > &Vec, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
Upwind the Fwd and Bwd states based on the velocity field given by Vec.
boost::shared_ptr< SegExp > SegExpSharedPtr
Base class for all multi-elemental spectral/hp expansions.
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
SpatialDomains::Geometry2DSharedPtr GetGeom2D() const
boost::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
boost::shared_ptr< SegGeom > SegGeomSharedPtr
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
void PeriodicEval(Array< OneD, NekDouble > &inarray1, Array< OneD, NekDouble > &inarray2, NekDouble h, int nmodes, Array< OneD, NekDouble > &outarray)
Evaluates the global spectral/hp expansion at some arbitray set of points.
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
int GetNumPoints() const
Return points order at which basis is defined.
boost::shared_ptr< Kernel > KernelSharedPtr
void SetCoeffPhysOffsets()
Definition of the total number of degrees of freedom and quadrature points and offsets to access data...
PointsManagerT & PointsManager(void)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
Defines a specification for a set of points.
boost::shared_ptr< GeometryVector > Composite
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
std::map< int, Composite > CompositeMap
boost::shared_ptr< Geometry2D > Geometry2DSharedPtr
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
void v_GetNormals(Array< OneD, Array< OneD, NekDouble > > &normals)
Populate normals with the normals of all expansions.
void Interp1D(const BasisKey &fbasis0, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, Array< OneD, NekDouble > &to)
this function interpolates a 1D function evaluated at the quadrature points of the basis fbasis0 to ...
This class is the abstraction of a one-dimensional multi-elemental expansions which is merely a colle...
boost::shared_ptr< Geometry1D > Geometry1DSharedPtr
void PostProcess(LibUtilities::KernelSharedPtr kernel, Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray, NekDouble h, int elmId=0)
Performs the post-processing on a specified element.
virtual void v_ReadGlobalOptimizationParameters()
void ReadGlobalOptimizationParameters()
void CreateCollections(Collections::ImplementationType ImpType=Collections::eNoImpType)
Construct collections of elements containing a single element type and polynomial order from the list...
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
int GetNumModes() const
Returns the order of the basis.
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
const StdRegions::StdExpansionVector &locexp);
void SetExpType(ExpansionType Type)
Returns the type of the expansion.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
boost::shared_ptr< StdExpansion > StdExpansionSharedPtr
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
boost::shared_ptr< MeshGraph > MeshGraphSharedPtr
boost::shared_ptr< Expansion1D > Expansion1DSharedPtr
int GetNedges() const
This function returns the number of edges of the expansion domain.
virtual ~ExpList1D()
Destructor.
Describes the specification for a Basis.
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
std::map< int, ExpansionShPtr > ExpansionMap
boost::shared_ptr< Expansion2D > Expansion2DSharedPtr