37 #include <boost/core/ignore_unused.hpp>
80 namespace MultiRegions
105 : m_expType(type), m_ncoeffs(0), m_npoints(0), m_physState(false),
121 : std::enable_shared_from_this<
ExpList>(in), m_expType(in.m_expType),
122 m_comm(in.m_comm), m_session(in.m_session), m_graph(in.m_graph),
123 m_ncoeffs(in.m_ncoeffs), m_npoints(in.m_npoints), m_physState(false),
124 m_exp(in.m_exp), m_collections(in.m_collections),
125 m_collectionsDoInit(in.m_collectionsDoInit),
126 m_coll_coeff_offset(in.m_coll_coeff_offset),
127 m_coll_phys_offset(in.m_coll_phys_offset),
128 m_coeff_offset(in.m_coeff_offset), m_phys_offset(in.m_phys_offset),
129 m_blockMat(in.m_blockMat), m_WaveSpace(false)
143 const bool DeclareCoeffPhysArrays,
145 : m_expType(in.m_expType), m_comm(in.m_comm), m_session(in.m_session),
146 m_graph(in.m_graph), m_physState(false),
151 for (
int i = 0; i < eIDs.size(); ++i)
153 (*m_exp).push_back((*(in.
m_exp))[eIDs[i]]);
185 const bool DeclareCoeffPhysArrays,
const std::string &var,
187 : m_comm(pSession->GetComm()), m_session(pSession), m_graph(graph),
195 graph->GetExpansionInfo(var);
228 const bool DeclareCoeffPhysArrays,
230 : m_comm(pSession->GetComm()), m_session(pSession), m_physState(false),
250 : m_expType(
e0D), m_ncoeffs(1), m_npoints(1), m_physState(false),
257 (*m_exp).push_back(
Point);
289 : m_comm(comm), m_session(pSession), m_graph(graph), m_physState(false),
294 boost::ignore_unused(variable, ImpType);
295 int i, j, id, elmtid = 0;
313 for (i = 0; i < bndCond.size(); ++i)
315 if (bndCond[i]->GetBoundaryConditionType() ==
318 for (j = 0; j < bndConstraint[i]->GetExpSize(); ++j)
321 std::dynamic_pointer_cast<LocalRegions::Expansion0D>(
322 bndConstraint[i]->
GetExp(j))))
326 PointGeom = exp0D->GetGeom()->GetVertex(0);
329 tracesDone.insert(PointGeom->GetVid());
331 else if ((exp1D = std::dynamic_pointer_cast<
333 bndConstraint[i]->
GetExp(j))))
338 exp1D->GetBasis(0)->GetBasisKey();
339 segGeom = exp1D->GetGeom1D();
343 tracesDone.insert(segGeom->GetGlobalID());
345 else if ((exp2D = std::dynamic_pointer_cast<
346 LocalRegions::Expansion2D>(
347 bndConstraint[i]->
GetExp(j))))
351 LibUtilities::BasisKey bkey0 =
352 exp2D->GetBasis(0)->GetBasisKey();
353 LibUtilities::BasisKey bkey1 =
354 exp2D->GetBasis(1)->GetBasisKey();
355 FaceGeom = exp2D->GetGeom2D();
358 if ((QuadGeom = std::dynamic_pointer_cast<
359 SpatialDomains::QuadGeom>(FaceGeom)))
362 LocalRegions::QuadExp>::AllocateSharedPtr(bkey0,
365 tracesDone.insert(QuadGeom->GetGlobalID());
368 else if ((TriGeom = std::dynamic_pointer_cast<
369 SpatialDomains::TriGeom>(FaceGeom)))
372 LocalRegions::TriExp>::AllocateSharedPtr(bkey0,
375 tracesDone.insert(TriGeom->GetGlobalID());
381 "proper face geometry failed");
385 exp->SetElmtId(elmtid++);
388 (*m_exp).push_back(exp);
393 map<int, pair<SpatialDomains::Geometry1DSharedPtr, LibUtilities::BasisKey>>
397 pair<LibUtilities::BasisKey, LibUtilities::BasisKey>>>
400 for (i = 0; i < locexp.size(); ++i)
402 if ((exp1D = std::dynamic_pointer_cast<LocalRegions::Expansion1D>(
407 for (j = 0; j < 2; ++j)
409 PointGeom = (exp1D->GetGeom1D())->GetVertex(j);
410 id = PointGeom->GetVid();
413 if (tracesDone.count(
id) != 0)
420 tracesDone.insert(
id);
421 exp->SetElmtId(elmtid++);
422 (*m_exp).push_back(exp);
425 else if ((exp2D = std::dynamic_pointer_cast<LocalRegions::Expansion2D>(
429 for (j = 0; j < locexp[i]->GetNtraces(); ++j)
431 segGeom = exp2D->GetGeom2D()->GetEdge(j);
432 id = segGeom->GetGlobalID();
434 if (tracesDone.count(
id) != 0)
439 auto it = edgeOrders.find(
id);
441 if (it == edgeOrders.end())
443 edgeOrders.insert(std::make_pair(
444 id, std::make_pair(segGeom,
445 locexp[i]->GetTraceBasisKey(j))));
449 LibUtilities::BasisKey edge =
450 locexp[i]->GetTraceBasisKey(j);
451 LibUtilities::BasisKey existing = it->second.second;
453 int np1 = edge.GetNumPoints();
454 int np2 = existing.GetNumPoints();
455 int nm1 = edge.GetNumModes();
456 int nm2 = existing.GetNumModes();
458 if (np2 >= np1 && nm2 >= nm1)
462 else if (np2 < np1 && nm2 < nm1)
464 it->second.second = edge;
469 "inappropriate number of points/modes (max"
470 "num of points is not set with max order)");
475 else if ((exp3D = dynamic_pointer_cast<LocalRegions::Expansion3D>(
479 for (j = 0; j < exp3D->GetNtraces(); ++j)
481 FaceGeom = exp3D->GetGeom3D()->GetFace(j);
482 id = FaceGeom->GetGlobalID();
484 if (tracesDone.count(
id) != 0)
488 auto it = faceOrders.find(
id);
490 if (it == faceOrders.end())
492 LibUtilities::BasisKey face_dir0 =
493 locexp[i]->GetTraceBasisKey(j, 0);
494 LibUtilities::BasisKey face_dir1 =
495 locexp[i]->GetTraceBasisKey(j, 1);
497 faceOrders.insert(std::make_pair(
499 std::make_pair(FaceGeom,
500 std::make_pair(face_dir0, face_dir1))));
504 LibUtilities::BasisKey face0 =
505 locexp[i]->GetTraceBasisKey(j, 0);
506 LibUtilities::BasisKey face1 =
507 locexp[i]->GetTraceBasisKey(j, 1);
508 LibUtilities::BasisKey existing0 = it->second.second.first;
509 LibUtilities::BasisKey existing1 = it->second.second.second;
511 int np11 = face0.GetNumPoints();
512 int np12 = face1.GetNumPoints();
513 int np21 = existing0.GetNumPoints();
514 int np22 = existing1.GetNumPoints();
515 int nm11 = face0.GetNumModes();
516 int nm12 = face1.GetNumModes();
517 int nm21 = existing0.GetNumModes();
518 int nm22 = existing1.GetNumModes();
520 if ((np22 >= np12 || np21 >= np11) &&
521 (nm22 >= nm12 || nm21 >= nm11))
525 else if ((np22 < np12 || np21 < np11) &&
526 (nm22 < nm12 || nm21 < nm11))
528 it->second.second.first = face0;
529 it->second.second.second = face1;
534 "inappropriate number of points/modes (max"
535 "num of points is not set with max order)");
542 int nproc =
m_comm->GetRowComm()->GetSize();
543 int tracepr =
m_comm->GetRowComm()->GetRank();
550 for (i = 0; i < locexp.size(); ++i)
552 tCnt += locexp[i]->GetNtraces();
557 Array<OneD, int> tracesCnt(nproc, 0);
558 tracesCnt[tracepr] = tCnt;
562 int totTraceCnt =
Vmath::Vsum(nproc, tracesCnt, 1);
563 Array<OneD, int> tTotOffsets(nproc, 0);
565 for (i = 1; i < nproc; ++i)
567 tTotOffsets[i] = tTotOffsets[i - 1] + tracesCnt[i - 1];
571 Array<OneD, int> TracesTotID(totTraceCnt, 0);
572 Array<OneD, int> TracesTotNm0(totTraceCnt, 0);
573 Array<OneD, int> TracesTotNm1(totTraceCnt, 0);
574 Array<OneD, int> TracesTotPnts0(totTraceCnt, 0);
575 Array<OneD, int> TracesTotPnts1(totTraceCnt, 0);
577 int cntr = tTotOffsets[tracepr];
579 for (i = 0; i < locexp.size(); ++i)
581 if ((exp2D = locexp[i]->as<LocalRegions::Expansion2D>()))
584 int nedges = locexp[i]->GetNtraces();
586 for (j = 0; j < nedges; ++j, ++cntr)
588 LibUtilities::BasisKey bkeyEdge =
589 locexp[i]->GetTraceBasisKey(j);
590 TracesTotID[cntr] = exp2D->GetGeom2D()->GetEid(j);
591 TracesTotNm0[cntr] = bkeyEdge.GetNumModes();
592 TracesTotPnts0[cntr] = bkeyEdge.GetNumPoints();
595 else if ((exp3D = locexp[i]->as<LocalRegions::Expansion3D>()))
597 int nfaces = locexp[i]->GetNtraces();
599 for (j = 0; j < nfaces; ++j, ++cntr)
601 LibUtilities::BasisKey face_dir0 =
602 locexp[i]->GetTraceBasisKey(j, 0);
603 LibUtilities::BasisKey face_dir1 =
604 locexp[i]->GetTraceBasisKey(j, 1);
606 TracesTotID[cntr] = exp3D->GetGeom3D()->GetFid(j);
607 TracesTotNm0[cntr] = face_dir0.GetNumModes();
608 TracesTotNm1[cntr] = face_dir1.GetNumModes();
609 TracesTotPnts0[cntr] = face_dir0.GetNumPoints();
610 TracesTotPnts1[cntr] = face_dir1.GetNumPoints();
617 m_comm->GetRowComm()->AllReduce(TracesTotPnts0,
625 if (edgeOrders.size())
627 for (i = 0; i < totTraceCnt; ++i)
629 auto it = edgeOrders.find(TracesTotID[i]);
631 if (it == edgeOrders.end())
636 LibUtilities::BasisKey existing = it->second.second;
637 LibUtilities::BasisKey edge(
638 existing.GetBasisType(), TracesTotNm0[i],
639 LibUtilities::PointsKey(TracesTotPnts0[i],
640 existing.GetPointsType()));
642 int np1 = edge.GetNumPoints();
643 int np2 = existing.GetNumPoints();
644 int nm1 = edge.GetNumModes();
645 int nm2 = existing.GetNumModes();
647 if (np2 >= np1 && nm2 >= nm1)
651 else if (np2 < np1 && nm2 < nm1)
653 it->second.second = edge;
658 "inappropriate number of points/modes (max "
659 "num of points is not set with max order)");
663 else if (faceOrders.size())
665 for (i = 0; i < totTraceCnt; ++i)
667 auto it = faceOrders.find(TracesTotID[i]);
669 if (it == faceOrders.end())
674 LibUtilities::BasisKey existing0 = it->second.second.first;
675 LibUtilities::BasisKey existing1 = it->second.second.second;
676 LibUtilities::BasisKey face0(
677 existing0.GetBasisType(), TracesTotNm0[i],
678 LibUtilities::PointsKey(TracesTotPnts0[i],
679 existing0.GetPointsType()));
680 LibUtilities::BasisKey face1(
681 existing1.GetBasisType(), TracesTotNm1[i],
682 LibUtilities::PointsKey(TracesTotPnts1[i],
683 existing1.GetPointsType()));
685 int np11 = face0.GetNumPoints();
686 int np12 = face1.GetNumPoints();
687 int np21 = existing0.GetNumPoints();
688 int np22 = existing1.GetNumPoints();
689 int nm11 = face0.GetNumModes();
690 int nm12 = face1.GetNumModes();
691 int nm21 = existing0.GetNumModes();
692 int nm22 = existing1.GetNumModes();
694 if ((np22 >= np12 || np21 >= np11) &&
695 (nm22 >= nm12 || nm21 >= nm11))
699 else if ((np22 < np12 || np21 < np11) &&
700 (nm22 < nm12 || nm21 < nm11))
702 it->second.second.first = face0;
703 it->second.second.second = face1;
708 "inappropriate number of points/modes (max "
709 "num of points is not set with max order)");
715 if (edgeOrders.size())
717 for (
auto &it : edgeOrders)
720 it.second.second, it.second.first);
721 exp->SetElmtId(elmtid++);
722 (*m_exp).push_back(exp);
727 for (
auto &it : faceOrders)
729 FaceGeom = it.second.first;
731 if ((QuadGeom = std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
735 it.second.second.first, it.second.second.second, QuadGeom);
738 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
742 it.second.second.first, it.second.second.second, TriGeom);
744 exp->SetElmtId(elmtid++);
745 (*m_exp).push_back(exp);
775 const bool DeclareCoeffPhysArrays,
const std::string variable,
777 : m_comm(pSession->GetComm()), m_session(pSession), m_graph(graph),
783 boost::ignore_unused(variable, ImpType);
784 int i, j, elmtid = 0;
799 for (i = 0; i < locexp.size(); ++i)
801 if ((exp1D = std::dynamic_pointer_cast<LocalRegions::Expansion1D>(
806 for (j = 0; j < 2; ++j)
808 PointGeom = (exp1D->GetGeom1D())->GetVertex(j);
812 exp->SetElmtId(elmtid++);
813 (*m_exp).push_back(exp);
816 else if ((exp2D = std::dynamic_pointer_cast<LocalRegions::Expansion2D>(
821 locexp[i]->GetBasis(0)->GetBasisKey();
823 for (j = 0; j < locexp[i]->GetNtraces(); ++j)
825 segGeom = exp2D->GetGeom2D()->GetEdge(j);
827 int dir = exp2D->GetGeom2D()->GetDir(j);
829 if (locexp[i]->GetNtraces() == 3)
832 locexp[i]->GetBasis(dir)->GetBasisKey();
846 locexp[i]->GetBasis(dir)->GetBasisKey(), segGeom);
849 exp->SetElmtId(elmtid++);
850 (*m_exp).push_back(exp);
853 else if ((exp3D = dynamic_pointer_cast<LocalRegions::Expansion3D>(
859 locexp[i]->GetBasis(0)->GetBasisKey();
861 locexp[i]->GetBasis(1)->GetBasisKey();
863 for (j = 0; j < exp3D->GetNtraces(); ++j)
865 FaceGeom = exp3D->GetGeom3D()->GetFace(j);
867 int dir0 = exp3D->GetGeom3D()->GetDir(j, 0);
868 int dir1 = exp3D->GetGeom3D()->GetDir(j, 1);
871 locexp[i]->GetBasis(dir0)->GetBasisKey();
873 locexp[i]->GetBasis(dir1)->GetBasisKey();
876 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
881 face_dir0, face_dir1, QuadGeom);
883 else if ((TriGeom = std::dynamic_pointer_cast<
895 nface_dir0, nface_dir1, TriGeom);
897 exp->SetElmtId(elmtid++);
898 (*m_exp).push_back(exp);
938 const bool DeclareCoeffPhysArrays,
const std::string variable,
939 bool SetToOneSpaceDimension,
942 : m_comm(comm), m_session(pSession), m_graph(graph), m_physState(false),
957 int meshdim = graph->GetMeshDimension();
961 graph->GetExpansionInfo(variable);
965 for (
auto &compIt : domain)
968 for (j = 0; j < compIt.second->m_geomVec.size(); ++j)
970 if ((PtGeom = std::dynamic_pointer_cast<SpatialDomains::PointGeom>(
971 compIt.second->m_geomVec[j])))
979 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
980 compIt.second->m_geomVec[j])))
989 auto expIt = expansions.find(SegGeom->GetGlobalID());
991 "Failed to find basis key");
992 bkey = expIt->second->m_basisKeyVector[0];
996 bkey = graph->GetEdgeBasisKey(SegGeom, variable);
999 if (SetToOneSpaceDimension)
1002 SegGeom->GenerateOneSpaceDimGeom();
1006 bkey, OneDSegmentGeom);
1017 std::dynamic_pointer_cast<SpatialDomains::TriGeom>(
1018 compIt.second->m_geomVec[j])))
1023 graph->GetFaceBasisKey(TriGeom, 0, variable);
1025 graph->GetFaceBasisKey(TriGeom, 1, variable);
1027 if (graph->GetExpansionInfo()
1029 ->second->m_basisKeyVector[0]
1045 TriBa, TriBb, TriGeom);
1048 else if ((QuadGeom =
1049 std::dynamic_pointer_cast<SpatialDomains::QuadGeom>(
1050 compIt.second->m_geomVec[j])))
1055 graph->GetFaceBasisKey(QuadGeom, 0, variable);
1057 graph->GetFaceBasisKey(QuadGeom, 1, variable);
1060 QuadBa, QuadBb, QuadGeom);
1065 "dynamic cast to a Geom (possibly 3D) failed");
1068 exp->SetElmtId(elmtid++);
1069 (*m_exp).push_back(exp);
1104 for (i = 0; i <
m_exp->size(); ++i)
1109 m_npoints += (*m_exp)[i]->GetTotPoints();
1113 if (DeclareCoeffPhysArrays)
1121 for (
int i = 0; i <
m_exp->size(); ++i)
1125 int loccoeffs = (*m_exp)[i]->GetNcoeffs();
1127 for (
int j = 0; j < loccoeffs; ++j)
1152 for (
auto &expIt : expmap)
1156 switch (expInfo->m_basisKeyVector.size())
1161 "Cannot mix expansion dimensions in one vector");
1165 std::dynamic_pointer_cast<SpatialDomains::SegGeom>(
1166 expInfo->m_geomShPtr)))
1178 "dynamic cast to a 1D Geom failed");
1185 "Cannot mix expansion dimensions in one vector");
1192 std::dynamic_pointer_cast<SpatialDomains ::TriGeom>(
1193 expInfo->m_geomShPtr)))
1214 else if ((QuadrilateralGeom = std::dynamic_pointer_cast<
1219 Ba, Bb, QuadrilateralGeom);
1224 "dynamic cast to a 2D Geom failed");
1231 "Cannot mix expansion dimensions in one vector");
1239 std::dynamic_pointer_cast<SpatialDomains::TetGeom>(
1240 expInfo->m_geomShPtr)))
1247 "LocalRegions::NodalTetExp is not implemented "
1257 else if ((PrismGeom = std::dynamic_pointer_cast<
1258 SpatialDomains ::PrismGeom>(
1259 expInfo->m_geomShPtr)))
1265 else if ((PyrGeom = std::dynamic_pointer_cast<
1271 Ba, Bb, Bc, PyrGeom);
1273 else if ((HexGeom = std::dynamic_pointer_cast<
1278 Ba, Bb, Bc, HexGeom);
1283 "dynamic cast to a Geom failed");
1289 "Dimension of basis key is greater than 3");
1293 exp->SetElmtId(
id++);
1296 (*m_exp).push_back(exp);
1326 int nrows = blockmat->GetRows();
1327 int ncols = blockmat->GetColumns();
1334 out = (*blockmat) * in;
1346 for (
int i = 0; i < (*m_exp).size(); ++i)
1348 (*m_exp)[i]->MultiplyByQuadratureMetric(inarray +
m_phys_offset[i],
1349 e_outarray = outarray +
1363 for (
int i = 0; i < (*m_exp).size(); ++i)
1365 (*m_exp)[i]->DivideByQuadratureMetric(inarray +
m_phys_offset[i],
1392 for (i = 0; i < (*m_exp).size(); ++i)
1394 (*m_exp)[i]->IProductWRTDerivBase(dir, inarray +
m_phys_offset[i],
1410 int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
1411 int nq = direction.size() / coordim;
1418 for (
int i = 0; i < (*m_exp).size(); ++i)
1420 npts_e = (*m_exp)[i]->GetTotPoints();
1423 for (
int k = 0; k < coordim; ++k)
1426 &locdir[k * npts_e], 1);
1429 (*m_exp)[i]->IProductWRTDirectionalDerivBase(
1454 ASSERTL1(inarray.size() >= dim,
"inarray is not of sufficient dimension");
1573 e_out_d0 = out_d0 + offset;
1574 e_out_d1 = out_d1 + offset;
1575 e_out_d2 = out_d2 + offset;
1578 inarray + offset, e_out_d0, e_out_d1,
1602 for (i = 0; i < (*m_exp).size(); ++i)
1605 (*m_exp)[i]->PhysDeriv_s(inarray +
m_phys_offset[i], e_out_ds);
1611 for (i = 0; i < (*m_exp).size(); i++)
1614 (*m_exp)[i]->PhysDeriv_n(inarray +
m_phys_offset[i], e_out_dn);
1631 int intdir = (int)edir;
1637 e_out_d = out_d + offset;
1640 inarray + offset, e_out_d);
1653 bool halfMode =
false;
1656 m_session->MatchSolverInfo(
"ModeType",
"HalfMode", halfMode,
false);
1708 ASSERTL0(0,
"Dimension not supported");
1719 int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
1720 int nq = direction.size() / coordim;
1726 for (
int i = 0; i < (*m_exp).size(); ++i)
1728 npts_e = (*m_exp)[i]->GetTotPoints();
1731 for (
int k = 0; k < coordim; ++k)
1734 &locdir[k * npts_e], 1);
1737 (*m_exp)[i]->PhysDirectionalDeriv(inarray +
m_phys_offset[i], locdir,
1749 for (
int i = 0; i < (*m_exp).size(); ++i)
1751 (*m_exp)[i]->ExponentialFilter(e_array = array +
m_phys_offset[i],
1752 alpha, exponent, cutoff);
1772 if (inarray.get() == outarray.get())
1775 out = (*InvMass) * in;
1780 out = (*InvMass) * in;
1819 for (i = 0; i < (*m_exp).size(); ++i)
1821 (*m_exp)[i]->FwdTransBndConstrained(inarray +
m_phys_offset[i],
1836 boost::ignore_unused(field);
1849 "Smoothing is currently not allowed unless you are using "
1850 "a nodal base for efficiency reasons. The implemented "
1851 "smoothing technique requires the mass matrix inversion "
1852 "which is trivial just for GLL_LAGRANGE_SEM and "
1853 "GAUSS_LAGRANGE_SEMexpansions.");
1881 map<int, int> elmt_id;
1886 for (i = 0; i < (*m_exp).size(); ++i)
1890 elmt_id[n_exp++] = i;
1896 n_exp = (*m_exp).size();
1897 for (i = 0; i < n_exp; ++i)
1911 for (i = 0; i < n_exp; ++i)
1913 nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
1914 ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1921 for (i = 0; i < n_exp; ++i)
1923 nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1924 ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetTotPoints();
1935 for (i = 0; i < n_exp; ++i)
1937 nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1938 ncols[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1946 for (i = 0; i < n_exp; ++i)
1948 nrows[i] = (*m_exp)[elmt_id.find(i)->second]->GetNcoeffs();
1950 (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
1957 for (i = 0; i < n_exp; ++i)
1960 (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
1962 (*m_exp)[elmt_id.find(i)->second]->NumDGBndryCoeffs();
1969 "Global Matrix creation not defined for this "
1982 for (i = cnt1 = 0; i < n_exp; ++i)
2001 loc_mat = std::dynamic_pointer_cast<LocalRegions::Expansion>(
2002 (*
m_exp)[elmt_id.find(i)->second])
2003 ->GetLocMatrix(matkey);
2004 BlkMatrix->SetBlock(i, i, loc_mat);
2021 return matrixIter->second;
2094 for (
int i = 0; i < (*m_exp).size(); ++i)
2112 (*m_exp)[i]->GeneralMatrixOp(
2129 int i, j, n, gid1, gid2, cntdim1, cntdim2;
2133 unsigned int glob_rows = 0;
2134 unsigned int glob_cols = 0;
2135 unsigned int loc_rows = 0;
2136 unsigned int loc_cols = 0;
2138 bool assembleFirstDim =
false;
2139 bool assembleSecondDim =
false;
2146 glob_cols = locToGloMap->GetNumGlobalCoeffs();
2148 assembleFirstDim =
false;
2149 assembleSecondDim =
true;
2154 glob_rows = locToGloMap->GetNumGlobalCoeffs();
2157 assembleFirstDim =
true;
2158 assembleSecondDim =
false;
2166 glob_rows = locToGloMap->GetNumGlobalCoeffs();
2167 glob_cols = locToGloMap->GetNumGlobalCoeffs();
2169 assembleFirstDim =
true;
2170 assembleSecondDim =
true;
2176 "Global Matrix creation not defined for this "
2188 for (n = cntdim1 = cntdim2 = 0; n < (*m_exp).size(); ++n)
2208 std::dynamic_pointer_cast<LocalRegions::Expansion>((*
m_exp)[n])
2209 ->GetLocMatrix(matkey);
2211 loc_rows = loc_mat->GetRows();
2212 loc_cols = loc_mat->GetColumns();
2214 for (i = 0; i < loc_rows; ++i)
2216 if (assembleFirstDim)
2218 gid1 = locToGloMap->GetLocalToGlobalMap(cntdim1 + i);
2219 sign1 = locToGloMap->GetLocalToGlobalSign(cntdim1 + i);
2227 for (j = 0; j < loc_cols; ++j)
2229 if (assembleSecondDim)
2231 gid2 = locToGloMap->GetLocalToGlobalMap(cntdim2 + j);
2232 sign2 = locToGloMap->GetLocalToGlobalSign(cntdim2 + j);
2241 coord = make_pair(gid1, gid2);
2242 if (spcoomat.count(coord) == 0)
2244 spcoomat[coord] = sign1 * sign2 * (*loc_mat)(i, j);
2248 spcoomat[coord] += sign1 * sign2 * (*loc_mat)(i, j);
2252 cntdim1 += loc_rows;
2253 cntdim2 += loc_cols;
2257 glob_cols, spcoomat);
2263 int i, j, n, gid1, gid2, loc_lda, eid;
2267 int totDofs = locToGloMap->GetNumGlobalCoeffs();
2268 int NumDirBCs = locToGloMap->GetNumGlobalDirBndCoeffs();
2270 unsigned int rows = totDofs - NumDirBCs;
2271 unsigned int cols = totDofs - NumDirBCs;
2275 int bwidth = locToGloMap->GetFullSystemBandWidth();
2287 if ((2 * (bwidth + 1)) < rows)
2291 rows, cols, zero, matStorage, bwidth, bwidth);
2297 rows, cols, zero, matStorage);
2311 for (n = 0; n < (*m_exp).size(); ++n)
2331 std::dynamic_pointer_cast<LocalRegions::Expansion>((*
m_exp)[n])
2332 ->GetLocMatrix(matkey);
2339 int rows = loc_mat->GetRows();
2340 int cols = loc_mat->GetColumns();
2341 const NekDouble *dat = loc_mat->GetRawPtr();
2344 Blas::Dscal(rows * cols, loc_mat->Scale(), new_mat->GetRawPtr(), 1);
2349 (*m_exp)[n]->AddRobinMassMatrix(
2350 rBC->m_robinID, rBC->m_robinPrimitiveCoeffs, new_mat);
2359 loc_lda = loc_mat->GetColumns();
2361 for (i = 0; i < loc_lda; ++i)
2365 sign1 = locToGloMap->GetLocalToGlobalSign(
m_coeff_offset[n] + i);
2368 for (j = 0; j < loc_lda; ++j)
2373 sign2 = locToGloMap->GetLocalToGlobalSign(
2381 if ((matStorage ==
eFULL) || (gid2 >= gid1))
2383 value = Gmat->GetValue(gid1, gid2) +
2384 sign1 * sign2 * (*loc_mat)(i, j);
2385 Gmat->SetValue(gid1, gid2, value);
2435 const map<int, RobinBCInfoSharedPtr> vRobinBCInfo =
GetRobinBCInfo();
2519 NekDouble tol,
bool returnNearestElmt,
int cachedId,
2524 return GetExpIndex(gloCoord, Lcoords, tol, returnNearestElmt, cachedId,
2530 bool returnNearestElmt,
int cachedId,
2543 for (
int i = (*m_exp).size() - 1; i >= 0; --i)
2554 if (cachedId >= 0 && cachedId < (*m_exp).size())
2557 if ((*
m_exp)[cachedId]->GetGeom()->MinMaxCheck(gloCoords) &&
2558 (*
m_exp)[cachedId]->GetGeom()->ContainsPoint(gloCoords, locCoords,
2563 else if (returnNearestElmt && (nearpt < nearpt_min))
2568 nearpt_min = nearpt;
2569 Vmath::Vcopy(locCoords.size(), locCoords, 1, savLocCoords, 1);
2573 NekDouble x = (gloCoords.size() > 0 ? gloCoords[0] : 0.0);
2574 NekDouble y = (gloCoords.size() > 1 ? gloCoords[1] : 0.0);
2575 NekDouble z = (gloCoords.size() > 2 ? gloCoords[2] : 0.0);
2582 std::vector<int> elmts =
m_graph->GetElementsContainingPoint(
p);
2585 for (
int i = 0; i < elmts.size(); ++i)
2592 if ((*
m_exp)[
id]->GetGeom()->ContainsPoint(gloCoords, locCoords, tol,
2597 else if (returnNearestElmt && (nearpt < nearpt_min))
2602 nearpt_min = nearpt;
2603 Vmath::Vcopy(locCoords.size(), locCoords, 1, savLocCoords, 1);
2609 if (returnNearestElmt && nearpt_min <= maxDistance)
2613 "Failed to find point within element to "
2615 boost::lexical_cast<std::string>(tol) +
" using local point (" +
2616 boost::lexical_cast<std::string>(locCoords[0]) +
"," +
2617 boost::lexical_cast<std::string>(locCoords[1]) +
"," +
2618 boost::lexical_cast<std::string>(locCoords[1]) +
2619 ") in element: " + boost::lexical_cast<std::string>(min_id);
2622 Vmath::Vcopy(locCoords.size(), savLocCoords, 1, locCoords, 1);
2639 ASSERTL0(dim == coords.size(),
"Invalid coordinate dimension.");
2644 ASSERTL0(elmtIdx > 0,
"Unable to find element containing point.");
2650 return (*
m_exp)[elmtIdx]->StdPhysEvaluate(xi, elmtPhys);
2679 for (
int i = 0; i <
m_exp->size(); ++i)
2681 (*m_exp)[i]->GetGeom()->Reset(
m_graph->GetCurvedEdges(),
2686 for (
int i = 0; i <
m_exp->size(); ++i)
2688 (*m_exp)[i]->Reset();
2704 int coordim =
GetExp(0)->GetCoordim();
2705 char vars[3] = {
'x',
'y',
'z'};
2716 outfile <<
"Variables = x";
2717 for (
int i = 1; i < coordim; ++i)
2719 outfile <<
", " << vars[i];
2724 outfile <<
", " << var;
2727 outfile << std::endl << std::endl;
2740 int nBases = (*m_exp)[0]->GetNumBases();
2745 if (expansion == -1)
2753 GetCoords(coords[0], coords[1], coords[2]);
2755 for (i = 0; i <
m_exp->size(); ++i)
2759 for (j = 0; j < nBases; ++j)
2761 numInt *= (*m_exp)[i]->GetNumPoints(j) - 1;
2764 numBlocks += numInt;
2769 nPoints = (*m_exp)[expansion]->GetTotPoints();
2775 (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
2778 for (j = 0; j < nBases; ++j)
2780 numBlocks *= (*m_exp)[expansion]->GetNumPoints(j) - 1;
2788 int nPlanes =
GetZIDs().size();
2789 NekDouble tmp = numBlocks * (nPlanes - 1.0) / nPlanes;
2790 numBlocks = (int)tmp;
2798 outfile <<
"Zone, N=" << nPoints <<
", E=" << numBlocks <<
", F=FEBlock";
2803 outfile <<
", ET=QUADRILATERAL" << std::endl;
2806 outfile <<
", ET=BRICK" << std::endl;
2814 for (j = 0; j < coordim; ++j)
2816 for (i = 0; i < nPoints; ++i)
2818 outfile << coords[j][i] <<
" ";
2819 if (i % 1000 == 0 && i)
2821 outfile << std::endl;
2824 outfile << std::endl;
2831 int nbase = (*m_exp)[0]->GetNumBases();
2834 std::shared_ptr<LocalRegions::ExpansionVector> exp =
m_exp;
2836 if (expansion != -1)
2838 exp = std::shared_ptr<LocalRegions::ExpansionVector>(
2840 (*exp)[0] = (*m_exp)[expansion];
2845 for (i = 0; i < (*exp).size(); ++i)
2847 const int np0 = (*exp)[i]->GetNumPoints(0);
2848 const int np1 = (*exp)[i]->GetNumPoints(1);
2850 for (j = 1; j < np1; ++j)
2852 for (k = 1; k < np0; ++k)
2854 outfile << cnt + (j - 1) * np0 + k <<
" ";
2855 outfile << cnt + (j - 1) * np0 + k + 1 <<
" ";
2856 outfile << cnt + j * np0 + k + 1 <<
" ";
2857 outfile << cnt + j * np0 + k << endl;
2864 else if (nbase == 3)
2866 for (i = 0; i < (*exp).size(); ++i)
2868 const int np0 = (*exp)[i]->GetNumPoints(0);
2869 const int np1 = (*exp)[i]->GetNumPoints(1);
2870 const int np2 = (*exp)[i]->GetNumPoints(2);
2871 const int np01 = np0 * np1;
2873 for (j = 1; j < np2; ++j)
2875 for (k = 1; k < np1; ++k)
2877 for (l = 1; l < np0; ++l)
2879 outfile << cnt + (j - 1) * np01 + (k - 1) * np0 + l
2881 outfile << cnt + (j - 1) * np01 + (k - 1) * np0 + l + 1
2883 outfile << cnt + (j - 1) * np01 + k * np0 + l + 1
2885 outfile << cnt + (j - 1) * np01 + k * np0 + l <<
" ";
2886 outfile << cnt + j * np01 + (k - 1) * np0 + l <<
" ";
2887 outfile << cnt + j * np01 + (k - 1) * np0 + l + 1
2889 outfile << cnt + j * np01 + k * np0 + l + 1 <<
" ";
2890 outfile << cnt + j * np01 + k * np0 + l << endl;
2894 cnt += np0 * np1 * np2;
2910 if (expansion == -1)
2918 for (
int i = 0; i < totpoints; ++i)
2920 outfile <<
m_phys[i] <<
" ";
2921 if (i % 1000 == 0 && i)
2923 outfile << std::endl;
2926 outfile << std::endl;
2930 int nPoints = (*m_exp)[expansion]->GetTotPoints();
2932 for (
int i = 0; i < nPoints; ++i)
2937 outfile << std::endl;
2943 outfile <<
"<?xml version=\"1.0\"?>" << endl;
2944 outfile <<
"<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" "
2945 <<
"byte_order=\"LittleEndian\">" << endl;
2946 outfile <<
" <UnstructuredGrid>" << endl;
2951 outfile <<
" </UnstructuredGrid>" << endl;
2952 outfile <<
"</VTKFile>" << endl;
2958 boost::ignore_unused(istrip);
2960 int nbase = (*m_exp)[expansion]->GetNumBases();
2961 int ntot = (*m_exp)[expansion]->GetTotPoints();
2965 for (i = 0; i < nbase; ++i)
2967 nquad[i] = (*m_exp)[expansion]->GetNumPoints(i);
2968 ntotminus *= (nquad[i] - 1);
2975 (*m_exp)[expansion]->GetCoords(coords[0], coords[1], coords[2]);
2977 outfile <<
" <Piece NumberOfPoints=\"" << ntot <<
"\" NumberOfCells=\""
2978 << ntotminus <<
"\">" << endl;
2979 outfile <<
" <Points>" << endl;
2980 outfile <<
" <DataArray type=\"Float64\" "
2981 <<
"NumberOfComponents=\"3\" format=\"ascii\">" << endl;
2983 for (i = 0; i < ntot; ++i)
2985 for (j = 0; j < 3; ++j)
2987 outfile << setprecision(8) << scientific << (float)coords[j][i]
2993 outfile <<
" </DataArray>" << endl;
2994 outfile <<
" </Points>" << endl;
2995 outfile <<
" <Cells>" << endl;
2996 outfile <<
" <DataArray type=\"Int32\" "
2997 <<
"Name=\"connectivity\" format=\"ascii\">" << endl;
3007 for (i = 0; i < nquad[0] - 1; ++i)
3009 outfile << i <<
" " << i + 1 << endl;
3017 for (i = 0; i < nquad[0] - 1; ++i)
3019 for (j = 0; j < nquad[1] - 1; ++j)
3021 outfile << j * nquad[0] + i <<
" " << j * nquad[0] + i + 1
3022 <<
" " << (j + 1) * nquad[0] + i + 1 <<
" "
3023 << (j + 1) * nquad[0] + i << endl;
3032 for (i = 0; i < nquad[0] - 1; ++i)
3034 for (j = 0; j < nquad[1] - 1; ++j)
3036 for (k = 0; k < nquad[2] - 1; ++k)
3039 << k * nquad[0] * nquad[1] + j * nquad[0] + i <<
" "
3040 << k * nquad[0] * nquad[1] + j * nquad[0] + i + 1
3042 << k * nquad[0] * nquad[1] + (j + 1) * nquad[0] +
3045 << k * nquad[0] * nquad[1] + (j + 1) * nquad[0] + i
3047 << (k + 1) * nquad[0] * nquad[1] + j * nquad[0] + i
3049 << (k + 1) * nquad[0] * nquad[1] + j * nquad[0] +
3052 << (k + 1) * nquad[0] * nquad[1] +
3053 (j + 1) * nquad[0] + i + 1
3055 << (k + 1) * nquad[0] * nquad[1] +
3056 (j + 1) * nquad[0] + i
3068 outfile <<
" </DataArray>" << endl;
3069 outfile <<
" <DataArray type=\"Int32\" "
3070 <<
"Name=\"offsets\" format=\"ascii\">" << endl;
3071 for (i = 0; i < ntotminus; ++i)
3073 outfile << i * ns + ns <<
" ";
3076 outfile <<
" </DataArray>" << endl;
3077 outfile <<
" <DataArray type=\"UInt8\" "
3078 <<
"Name=\"types\" format=\"ascii\">" << endl;
3079 for (i = 0; i < ntotminus; ++i)
3084 outfile <<
" </DataArray>" << endl;
3085 outfile <<
" </Cells>" << endl;
3086 outfile <<
" <PointData>" << endl;
3091 boost::ignore_unused(expansion);
3092 outfile <<
" </PointData>" << endl;
3093 outfile <<
" </Piece>" << endl;
3100 int nq = (*m_exp)[expansion]->GetTotPoints();
3103 outfile <<
" <DataArray type=\"Float64\" Name=\"" << var <<
"\">"
3107 for (i = 0; i < nq; ++i)
3113 outfile <<
" </DataArray>" << endl;
3147 err = max(err,
abs(inarray[i] - soln[i]));
3180 for (i = 0; i < (*m_exp).size(); ++i)
3183 err += errl2 * errl2;
3188 for (i = 0; i < (*m_exp).size(); ++i)
3192 err += errl2 * errl2;
3222 for (i = 0; i < (*m_exp).size(); ++i)
3238 for (i = 0; i < (*m_exp).size(); ++i)
3241 for (j = 0; j < inarray.size(); ++j)
3245 flux += (*m_exp)[i]->VectorFlux(tmp);
3254 "This method is not defined or valid for this class type");
3262 "This method is not defined or valid for this class type");
3270 "This method is not defined or valid for this class type");
3277 boost::ignore_unused(lhom);
3279 "This method is not defined or valid for this class type");
3285 "This method is not defined or valid for this class type");
3293 "This method is not defined or valid for this class type");
3301 "This method is not defined or valid for this class type");
3306 const std::string &varName,
3307 const std::shared_ptr<ExpList> locExpList)
3309 string varString = fileName.substr(0, fileName.find_last_of(
"."));
3310 int j, k, len = varString.length();
3311 varString = varString.substr(len - 1, len);
3313 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
3314 std::vector<std::vector<NekDouble>> FieldData;
3319 ft, comm,
m_session->GetSharedFilesystem());
3321 f->Import(fileName, FieldDef, FieldData);
3324 for (j = 0; j < FieldDef.size(); ++j)
3326 for (k = 0; k < FieldDef[j]->m_fields.size(); ++k)
3328 if (FieldDef[j]->m_fields[k] == varName)
3331 locExpList->ExtractDataToCoeffs(FieldDef[j], FieldData[j],
3332 FieldDef[j]->m_fields[k],
3333 locExpList->UpdateCoeffs());
3339 ASSERTL0(found,
"Could not find variable '" + varName +
3340 "' in file boundary condition " + fileName);
3341 locExpList->BwdTrans(locExpList->GetCoeffs(), locExpList->UpdatePhys());
3367 for (i = 0; i < (*m_exp).size(); ++i)
3371 err += errh1 * errh1;
3380 std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef,
3382 std::vector<NekDouble> &HomoLen,
bool homoStrips,
3383 std::vector<unsigned int> &HomoSIDs, std::vector<unsigned int> &HomoZIDs,
3384 std::vector<unsigned int> &HomoYIDs)
3391 ASSERTL1(NumHomoDir == HomoBasis.size(),
3392 "Homogeneous basis is not the same length as NumHomoDir");
3393 ASSERTL1(NumHomoDir == HomoLen.size(),
3394 "Homogeneous length vector is not the same length as NumHomDir");
3413 for (s = startenum; s <= endenum; ++s)
3415 std::vector<unsigned int> elementIDs;
3416 std::vector<LibUtilities::BasisType> basis;
3417 std::vector<unsigned int> numModes;
3418 std::vector<std::string> fields;
3421 bool UniOrder =
true;
3426 for (
int i = 0; i < (*m_exp).size(); ++i)
3428 if ((*
m_exp)[i]->GetGeom()->GetShapeType() == shape)
3430 elementIDs.push_back((*
m_exp)[i]->GetGeom()->GetGlobalID());
3433 for (
int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
3436 (*
m_exp)[i]->GetBasis(j)->GetBasisType());
3438 (*
m_exp)[i]->GetBasis(j)->GetNumModes());
3442 for (n = 0; n < NumHomoDir; ++n)
3444 basis.push_back(HomoBasis[n]->GetBasisType());
3445 numModes.push_back(HomoBasis[n]->GetNumModes());
3453 (*
m_exp)[i]->GetBasis(0)->GetBasisType() == basis[0],
3454 "Routine is not set up for multiple bases definitions");
3456 for (
int j = 0; j < (*m_exp)[i]->GetNumBases(); ++j)
3459 (*
m_exp)[i]->GetBasis(j)->GetNumModes());
3461 (*
m_exp)[i]->GetBasis(j)->GetNumModes())
3467 for (n = 0; n < NumHomoDir; ++n)
3469 numModes.push_back(HomoBasis[n]->GetNumModes());
3475 if (elementIDs.size() > 0)
3480 numModes, fields, NumHomoDir, HomoLen,
3481 homoStrips, HomoSIDs, HomoZIDs, HomoYIDs);
3482 fielddef.push_back(fdef);
3490 std::vector<LibUtilities::FieldDefinitionsSharedPtr>
ExpList::
3493 std::vector<LibUtilities::FieldDefinitionsSharedPtr> returnval;
3499 std::vector<LibUtilities::FieldDefinitionsSharedPtr> &fielddef)
3508 std::vector<NekDouble> &fielddata)
3522 map<int, int> ElmtID_to_ExpID;
3523 for (i = 0; i < (*m_exp).size(); ++i)
3525 ElmtID_to_ExpID[(*m_exp)[i]->GetGeom()->GetGlobalID()] = i;
3528 for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
3530 int eid = ElmtID_to_ExpID[fielddef->m_elementIDs[i]];
3531 int datalen = (*m_exp)[eid]->GetNcoeffs();
3540 std::vector<NekDouble> &fielddata, std::string &field,
3547 const std::shared_ptr<ExpList> &fromExpList,
3564 std::vector<NekDouble> &fielddata, std::string &field,
3569 int modes_offset = 0;
3570 int datalen = fielddata.size() / fielddef->m_fields.size();
3573 for (i = 0; i < fielddef->m_fields.size(); ++i)
3575 if (fielddef->m_fields[i] == field)
3582 ASSERTL0(i != fielddef->m_fields.size(),
3583 "Field (" + field +
") not found in file.");
3590 for (i = (*m_exp).size() - 1; i >= 0; --i)
3596 for (i = 0; i < fielddef->m_elementIDs.size(); ++i)
3600 if (fielddef->m_uniOrder ==
true)
3606 fielddef->m_shapeType, fielddef->m_numModes, modes_offset);
3608 const int elmtId = fielddef->m_elementIDs[i];
3614 modes_offset += (*m_exp)[0]->GetNumBases();
3618 expId = eIt->second;
3620 bool sameBasis =
true;
3621 for (
int j = 0; j < fielddef->m_basis.size(); ++j)
3623 if (fielddef->m_basis[j] != (*
m_exp)[expId]->GetBasisType(j))
3637 (*m_exp)[expId]->ExtractDataToCoeffs(
3638 &fielddata[offset], fielddef->m_numModes, modes_offset,
3643 modes_offset += (*m_exp)[0]->GetNumBases();
3650 const std::shared_ptr<ExpList> &fromExpList,
3657 for (i = 0; i < (*m_exp).size(); ++i)
3659 std::vector<unsigned int> nummodes;
3660 vector<LibUtilities::BasisType> basisTypes;
3661 for (
int j = 0; j < fromExpList->GetExp(i)->GetNumBases(); ++j)
3663 nummodes.push_back(fromExpList->GetExp(i)->GetBasisNumModes(j));
3664 basisTypes.push_back(fromExpList->GetExp(i)->GetBasisType(j));
3667 (*m_exp)[i]->ExtractDataToCoeffs(&fromCoeffs[offset], nummodes, 0,
3671 offset += fromExpList->GetExp(i)->GetNcoeffs();
3681 size_t nTracePts = weightAver.size();
3683 for (
int i = 0; i < nTracePts; ++i)
3685 weightAver[i] = 0.5;
3686 weightJump[i] = 1.0;
3698 int nq = outarray[0].size() / MFdim;
3701 int coordim = (*m_exp)[0]->GetGeom()->GetCoordim();
3705 for (
int i = 0; i <
m_exp->size(); ++i)
3707 npts = (*m_exp)[i]->GetTotPoints();
3709 for (
int j = 0; j < MFdim * coordim; ++j)
3715 (*m_exp)[i]->GetMetricInfo()->GetMovingFrames(
3716 (*
m_exp)[i]->GetPointsKeys(), MMFdir, CircCentre, MFloc);
3719 for (
int j = 0; j < MFdim; ++j)
3721 for (
int k = 0; k < coordim; ++k)
3745 for (
int i = 0; i < (*m_exp).size(); ++i)
3747 npoints_e = (*m_exp)[i]->GetTotPoints();
3768 "This method is not defined or valid for this class type");
3775 boost::ignore_unused(i);
3777 "This method is not defined or valid for this class type");
3778 static std::shared_ptr<ExpList> result;
3801 int i, j, k, e_npoints, offset;
3809 "Input vector does not have sufficient dimensions to "
3813 for (i = 0; i <
m_exp->size(); ++i)
3816 e_npoints = (*m_exp)[i]->GetNumPoints(0);
3817 normals = (*m_exp)[i]->GetPhysNormals();
3823 for (j = 0; j < e_npoints; ++j)
3827 for (k = 0; k < coordim; ++k)
3829 Vn += Vec[k][offset + j] * normals[k * e_npoints + j];
3835 Upwind[offset + j] = Fwd[offset + j];
3839 Upwind[offset + j] = Bwd[offset + j];
3847 "This method is not defined or valid for this class type");
3893 "This method is not defined or valid for this class type");
3894 static std::shared_ptr<ExpList> returnVal;
3901 "This method is not defined or valid for this class type");
3902 static std::shared_ptr<AssemblyMapDG> result;
3908 return GetTraceMap()->GetBndCondIDToGlobalTraceID();
3914 "This method is not defined or valid for this class type");
3915 static std::vector<bool> result;
3932 for (
int i = 0; i < nquad2; ++i)
3934 for (
int j = 0; j < nquad1; ++j)
3936 out[i * nquad1 + j] = in[j * nquad2 + i];
3942 for (
int i = 0; i < nquad2; ++i)
3944 for (
int j = 0; j < nquad1; ++j)
3946 out[i * nquad1 + j] = in[i * nquad1 + j];
3957 for (
int i = 0; i < nquad2; ++i)
3959 for (
int j = 0; j < nquad1 / 2; ++j)
3961 swap(out[i * nquad1 + j], out[i * nquad1 + nquad1 - j - 1]);
3972 for (
int j = 0; j < nquad1; ++j)
3974 for (
int i = 0; i < nquad2 / 2; ++i)
3976 swap(out[i * nquad1 + j], out[(nquad2 - i - 1) * nquad1 + j]);
3991 int i, j, k, e_npoints, offset;
3997 ASSERTL1(normals.size() >= coordim,
3998 "Output vector does not have sufficient dimensions to "
4006 for (i = 0; i <
m_exp->size(); ++i)
4011 loc_exp->GetLeftAdjacentElementExp();
4015 locnormals = loc_elmt->GetTraceNormal(
4016 loc_exp->GetLeftAdjacentElementTrace());
4022 for (j = 0; j < e_npoints; ++j)
4026 for (k = 0; k < coordim; ++k)
4028 normals[k][offset] = locnormals[k][0];
4040 for (i = 0; i <
m_exp->size(); ++i)
4045 loc_exp->GetLeftAdjacentElementExp();
4047 int edgeNumber = loc_exp->GetLeftAdjacentElementTrace();
4050 e_npoints = (*m_exp)[i]->GetNumPoints(0);
4052 locnormals = loc_elmt->GetTraceNormal(edgeNumber);
4053 int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
4054 int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
4056 if (e_nmodes != loc_nmodes)
4058 if (loc_exp->GetRightAdjacentElementTrace() >= 0)
4061 loc_exp->GetRightAdjacentElementExp();
4064 loc_exp->GetRightAdjacentElementTrace();
4068 locnormals = loc_elmt->GetTraceNormal(EdgeNumber);
4073 for (j = 0; j < e_npoints; ++j)
4078 for (k = 0; k < coordim; ++k)
4080 normals[k][offset + j] = -locnormals[k][j];
4089 for (
int p = 0;
p < coordim; ++
p)
4093 loc_exp->GetBasis(0)->GetPointsKey();
4095 loc_elmt->GetBasis(0)->GetPointsKey();
4103 for (j = 0; j < e_npoints; ++j)
4107 for (k = 0; k < coordim; ++k)
4109 normals[k][offset + j] = normal[k][j];
4120 for (j = 0; j < e_npoints; ++j)
4124 for (k = 0; k < coordim; ++k)
4126 normals[k][offset + j] = locnormals[k][j];
4138 for (i = 0; i <
m_exp->size(); ++i)
4142 traceExp->GetLeftAdjacentElementExp();
4145 int faceNum = traceExp->GetLeftAdjacentElementTrace();
4149 exp3D->GetTraceNormal(faceNum);
4155 int fromid0, fromid1;
4169 exp3D->GetTraceBasisKey(faceNum, fromid0);
4171 exp3D->GetTraceBasisKey(faceNum, fromid1);
4173 traceExp->GetBasis(0)->GetBasisKey();
4175 traceExp->GetBasis(1)->GetBasisKey();
4181 exp3D->ReOrientTracePhysMap(orient, faceids, faceNq0, faceNq1);
4185 for (j = 0; j < coordim; ++j)
4187 Vmath::Scatr(faceNq0 * faceNq1, locNormals[j], faceids,
4193 traceBasis1.
GetPointsKey(), tmp = normals[j] + offset);
4201 "This method is not defined or valid for this class type");
4216 lengLR[0] = lengthsFwd;
4217 lengLR[1] = lengthsBwd;
4221 int e_npoints0 = -1;
4224 for (
int i = 0; i <
m_exp->size(); ++i)
4226 loc_exp = (*m_exp)[i];
4229 int e_nmodes = loc_exp->GetBasis(0)->GetNumModes();
4230 e_npoints = (*m_exp)[i]->GetNumPoints(0);
4231 if (e_npoints0 < e_npoints)
4233 for (
int nlr = 0; nlr < 2; nlr++)
4237 e_npoints0 = e_npoints;
4240 LRelmts[0] = loc_exp->GetLeftAdjacentElementExp();
4241 LRelmts[1] = loc_exp->GetRightAdjacentElementExp();
4243 LRbndnumbs[0] = loc_exp->GetLeftAdjacentElementTrace();
4244 LRbndnumbs[1] = loc_exp->GetRightAdjacentElementTrace();
4245 for (
int nlr = 0; nlr < 2; ++nlr)
4248 lengAdd[nlr] = lengintp[nlr];
4249 int bndNumber = LRbndnumbs[nlr];
4250 loc_elmt = LRelmts[nlr];
4253 locLeng = loc_elmt->GetElmtBndNormDirElmtLen(bndNumber);
4254 lengAdd[nlr] = locLeng;
4256 int loc_nmodes = loc_elmt->GetBasis(0)->GetNumModes();
4257 if (e_nmodes != loc_nmodes)
4261 loc_exp->GetBasis(0)->GetPointsKey();
4263 loc_elmt->GetBasis(0)->GetPointsKey();
4266 lengAdd[nlr] = lengintp[nlr];
4270 for (
int j = 0; j < e_npoints; ++j)
4272 lengLR[nlr][offset + j] = lengAdd[nlr][j];
4279 for (
int i = 0; i <
m_exp->size(); ++i)
4281 loc_exp = (*m_exp)[i];
4285 loc_exp->GetBasis(0)->GetBasisKey();
4287 loc_exp->GetBasis(1)->GetBasisKey();
4290 e_npoints = TraceNq0 * TraceNq1;
4291 if (e_npoints0 < e_npoints)
4293 for (
int nlr = 0; nlr < 2; nlr++)
4297 e_npoints0 = e_npoints;
4300 LRelmts[0] = loc_exp->GetLeftAdjacentElementExp();
4301 LRelmts[1] = loc_exp->GetRightAdjacentElementExp();
4303 LRbndnumbs[0] = loc_exp->GetLeftAdjacentElementTrace();
4304 LRbndnumbs[1] = loc_exp->GetRightAdjacentElementTrace();
4305 for (
int nlr = 0; nlr < 2; ++nlr)
4308 int bndNumber = LRbndnumbs[nlr];
4309 loc_elmt = LRelmts[nlr];
4312 locLeng = loc_elmt->GetElmtBndNormDirElmtLen(bndNumber);
4316 loc_elmt->GetTraceOrient(bndNumber);
4318 int fromid0, fromid1;
4331 loc_elmt->GetTraceBasisKey(bndNumber, fromid0);
4333 loc_elmt->GetTraceBasisKey(bndNumber, fromid1);
4338 AlignFace(orient, faceNq0, faceNq1, locLeng, alignedLeng);
4345 for (
int j = 0; j < e_npoints; ++j)
4347 lengLR[nlr][offset + j] = lengintp[nlr][j];
4358 boost::ignore_unused(Fx, Fy, outarray);
4360 "This method is not defined or valid for this class type");
4366 boost::ignore_unused(Fn, outarray);
4368 "This method is not defined or valid for this class type");
4375 boost::ignore_unused(Fwd, Bwd, outarray);
4377 "This method is not defined or valid for this class type");
4383 boost::ignore_unused(Fwd, Bwd);
4385 "This method is not defined or valid for this class type");
4391 bool PutFwdInBwdOnBCs,
bool DoExchange)
4393 boost::ignore_unused(field, Fwd, Bwd, FillBnd, PutFwdInBwdOnBCs,
4396 "This method is not defined or valid for this class type");
4401 bool PutFwdInBwdOnBCs)
4403 boost::ignore_unused(Fwd, Bwd, PutFwdInBwdOnBCs);
4405 "This method is not defined or valid for this class type");
4412 boost::ignore_unused(field, Fwd, Bwd);
4414 "v_AddTraceQuadPhysToField is not defined for this class type");
4421 boost::ignore_unused(field, Fwd, Bwd);
4423 "v_AddTraceQuadPhysToOffDiag is not defined for this class");
4431 boost::ignore_unused(Fwd, Bwd, locTraceFwd, locTraceBwd);
4433 "v_GetLocTraceFromTracePts is not defined for this class");
4439 "v_GetBndCondBwdWeight is not defined for this class type");
4446 boost::ignore_unused(index, value);
4448 "v_setBndCondBwdWeight is not defined for this class type");
4454 "This method is not defined or valid for this class type");
4455 static vector<bool> tmp;
4461 boost::ignore_unused(outarray);
4463 "This method is not defined or valid for this class type");
4469 boost::ignore_unused(inarray, outarray);
4471 "This method is not defined or valid for this class type");
4478 boost::ignore_unused(inarray, outarray);
4480 "This method is not defined or valid for this class type");
4489 const bool PhysSpaceForcing)
4491 boost::ignore_unused(inarray, outarray, factors, varcoeff, varfactors,
4492 dirForcing, PhysSpaceForcing);
4502 boost::ignore_unused(velocity, inarray, outarray, lambda, dirForcing);
4504 "This method is not defined or valid for this class type");
4513 boost::ignore_unused(velocity, inarray, outarray, lambda, dirForcing);
4515 "This method is not defined or valid for this class type");
4520 bool Shuff,
bool UnShuff)
4522 boost::ignore_unused(inarray, outarray, Shuff, UnShuff);
4524 "This method is not defined or valid for this class type");
4529 bool Shuff,
bool UnShuff)
4531 boost::ignore_unused(inarray, outarray, Shuff, UnShuff);
4533 "This method is not defined or valid for this class type");
4540 boost::ignore_unused(inarray1, inarray2, outarray);
4542 "This method is not defined or valid for this class type");
4550 boost::ignore_unused(inarray1, inarray2, outarray);
4552 "This method is not defined or valid for this class type");
4558 boost::ignore_unused(BndVals, TotField, BndID);
4560 "This method is not defined or valid for this class type");
4568 boost::ignore_unused(V1, V2, outarray, BndID);
4570 "This method is not defined or valid for this class type");
4583 (*m_exp)[i]->NormVectorIProductWRTBase(
4593 (*m_exp)[i]->NormVectorIProductWRTBase(
4603 (*m_exp)[i]->NormVectorIProductWRTBase(
4618 boost::ignore_unused(outarray);
4620 "This method is not defined or valid for this class type");
4628 "This method is not defined or valid for this class type");
4635 boost::ignore_unused(nreg);
4637 "This method is not defined or valid for this class type");
4642 boost::ignore_unused(useComm);
4644 "This method is not defined or valid for this class type");
4650 boost::ignore_unused(inarray, outarray, useComm);
4652 "This method is not defined or valid for this class type");
4658 "This method is not defined or valid for this class type");
4664 boost::ignore_unused(inarray, outarray);
4666 "This method is not defined or valid for this class type");
4737 for (i = 0; i < (*m_exp).size(); ++i)
4740 (*m_exp)[i]->GetCoords(e_coord_0);
4744 ASSERTL0(coord_1.size() != 0,
"output coord_1 is not defined");
4746 for (i = 0; i < (*m_exp).size(); ++i)
4750 (*m_exp)[i]->GetCoords(e_coord_0, e_coord_1);
4754 ASSERTL0(coord_1.size() != 0,
"output coord_1 is not defined");
4755 ASSERTL0(coord_2.size() != 0,
"output coord_2 is not defined");
4757 for (i = 0; i < (*m_exp).size(); ++i)
4762 (*m_exp)[i]->GetCoords(e_coord_0, e_coord_1, e_coord_2);
4772 (*m_exp)[eid]->GetCoords(xc0, xc1, xc2);
4782 for (
int i = 0; i <
m_exp->size(); ++i)
4784 for (
int j = 0; j < (*m_exp)[i]->GetNtraces(); ++j)
4786 (*m_exp)[i]->ComputeTraceNormal(j);
4794 const bool DeclareCoeffPhysArrays)
4796 boost::ignore_unused(i, result, DeclareCoeffPhysArrays);
4798 "This method is not defined or valid for this class type");
4819 for (cnt = n = 0; n < i; ++n)
4830 elmt =
GetExp(ElmtID[cnt + n]);
4831 elmt->GetTracePhysVals(
4833 tmp1 = element + offsetElmt, tmp2 = boundary + offsetBnd);
4835 offsetElmt += elmt->GetTotPoints();
4851 for (cnt = n = 0; n < i; ++n)
4860 npoints +=
GetExp(ElmtID[cnt + n])->GetTotPoints();
4871 nq =
GetExp(ElmtID[cnt + n])->GetTotPoints();
4873 Vmath::Vcopy(nq, &phys[offsetPhys], 1, &bndElmt[offsetElmt], 1);
4896 for (cnt = n = 0; n < i; ++n)
4908 elmt =
GetExp(ElmtID[cnt + n]);
4909 elmt->GetTracePhysVals(EdgeID[cnt + n],
4911 phys + offsetPhys, tmp1 = bnd + offsetBnd);
4930 for (j = 0; j < coordim; ++j)
4937 for (cnt = n = 0; n < i; ++n)
4948 elmt =
GetExp(ElmtID[cnt + n]);
4950 elmt->GetTraceNormal(EdgeID[cnt + n]);
4952 for (j = 0; j < coordim; ++j)
4954 Vmath::Vcopy(nq, normalsElmt[j], 1, tmp = normals[j] + offset, 1);
4964 boost::ignore_unused(ElmtID, EdgeID);
4966 "This method is not defined or valid for this class type");
4972 boost::ignore_unused(weightave, weightjmp);
4979 boost::ignore_unused(Fwd, Bwd);
4989 "This method is not defined or valid for this class type");
5000 "This method is not defined or valid for this class type");
5008 const std::string varName,
5012 boost::ignore_unused(time, varName, x2_in, x3_in);
5014 "This method is not defined or valid for this class type");
5022 "This method is not defined or valid for this class type");
5023 static map<int, RobinBCInfoSharedPtr> result;
5033 boost::ignore_unused(periodicVerts, periodicEdges, periodicFaces);
5035 "This method is not defined or valid for this class type");
5040 unsigned int regionId,
const std::string &variable)
5042 auto collectionIter = collection.find(regionId);
5043 ASSERTL1(collectionIter != collection.end(),
5044 "Unable to locate collection " +
5045 boost::lexical_cast<string>(regionId));
5048 (*collectionIter).second;
5049 auto conditionMapIter = bndCondMap->find(variable);
5050 ASSERTL1(conditionMapIter != bndCondMap->end(),
5051 "Unable to locate condition map.");
5054 (*conditionMapIter).second;
5056 return boundaryCondition;
5061 boost::ignore_unused(n);
5063 "This method is not defined or valid for this class type");
5074 vector<std::pair<LocalRegions::ExpansionSharedPtr, int>>>
5100 if (
m_graph->GetMeshDimension() == (*
m_exp)[0]->GetShapeDimension())
5106 bool verbose = (
m_session->DefinesCmdLineArgument(
"verbose")) &&
5110 : 2 *
m_exp->size());
5118 for (
int i = 0; i <
m_exp->size(); ++i)
5120 collections[(*m_exp)[i]->DetShapeType()].push_back(
5121 std::pair<LocalRegions::ExpansionSharedPtr, int>((*
m_exp)[i], i));
5124 for (
auto &it : collections)
5130 vector<StdRegions::StdExpansionSharedPtr> collExp;
5141 if (it.second.size() == 1)
5143 collExp.push_back(it.second[0].first);
5150 if (collExp.size() > collsize)
5154 collsize = collExp.size();
5164 collExp.push_back(it.second[0].first);
5165 int prevnCoeff = it.second[0].first->GetNcoeffs();
5166 int prevnPhys = it.second[0].first->GetTotPoints();
5168 it.second[0].first->GetMetricInfo()->GetGtype() ==
5172 for (
int i = 1; i < it.second.size(); ++i)
5174 int nCoeffs = it.second[i].first->GetNcoeffs();
5175 int nPhys = it.second[i].first->GetTotPoints();
5177 it.second[i].first->GetMetricInfo()->GetGtype() ==
5185 if (prevCoeffOffset + nCoeffs != coeffOffset ||
5186 prevnCoeff != nCoeffs ||
5187 prevPhysOffset + nPhys != physOffset ||
5188 prevDeformed != Deformed || prevnPhys != nPhys ||
5197 if (collExp.size() > collsize)
5201 collsize = collExp.size();
5213 collExp.push_back(it.second[i].first);
5218 collExp.push_back(it.second[i].first);
5223 if (i == it.second.size() - 1)
5229 if (collExp.size() > collsize)
5233 collsize = collExp.size();
5244 prevCoeffOffset = coeffOffset;
5245 prevPhysOffset = physOffset;
5246 prevDeformed = Deformed;
5247 prevnCoeff = nCoeffs;
5283 int pt0 = (*m_exp)[i]->GetNumPoints(0);
5284 int pt1 = (*m_exp)[i]->GetNumPoints(1);
5285 int npt0 = (int)pt0 * scale;
5286 int npt1 = (int)pt1 * scale;
5289 npt0, (*
m_exp)[i]->GetPointsType(0));
5291 npt1, (*
m_exp)[i]->GetPointsType(1));
5295 (*
m_exp)[i]->GetBasis(1)->GetPointsKey(),
5296 &inarray[cnt], newPointsKey0,
5297 newPointsKey1, &outarray[cnt1]);
5300 cnt1 += npt0 * npt1;
5309 int pt0 = (*m_exp)[i]->GetNumPoints(0);
5310 int pt1 = (*m_exp)[i]->GetNumPoints(1);
5311 int pt2 = (*m_exp)[i]->GetNumPoints(2);
5312 int npt0 = (int)pt0 * scale;
5313 int npt1 = (int)pt1 * scale;
5314 int npt2 = (int)pt2 * scale;
5317 npt0, (*
m_exp)[i]->GetPointsType(0));
5319 npt1, (*
m_exp)[i]->GetPointsType(1));
5321 npt2, (*
m_exp)[i]->GetPointsType(2));
5325 (*
m_exp)[i]->GetBasis(1)->GetPointsKey(),
5326 (*
m_exp)[i]->GetBasis(2)->GetPointsKey(),
5327 &inarray[cnt], newPointsKey0,
5328 newPointsKey1, newPointsKey2,
5331 cnt += pt0 * pt1 * pt2;
5332 cnt1 += npt0 * npt1 * npt2;
5349 boost::ignore_unused(FwdFlux, BwdFlux, outarray);
5357 int nTotElmt = (*m_exp).size();
5358 int nElmtPnt = (*m_exp)[0]->GetTotPoints();
5359 int nElmtCoef = (*m_exp)[0]->GetNcoeffs();
5364 for (
int nelmt = 0; nelmt < nTotElmt; nelmt++)
5366 nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
5367 nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
5369 if (tmpPhys.size() != nElmtPnt || tmpCoef.size() != nElmtCoef)
5375 for (
int ncl = 0; ncl < nElmtPnt; ncl++)
5377 tmpPhys[ncl] = inarray[nelmt][ncl];
5379 (*m_exp)[nelmt]->IProductWRTDerivBase(nDirctn, tmpPhys, tmpCoef);
5381 for (
int nrw = 0; nrw < nElmtCoef; nrw++)
5383 (*mtxPerVar[nelmt])(nrw, ncl) = tmpCoef[nrw];
5394 int nTotElmt = (*m_exp).size();
5396 int nspacedim =
m_graph->GetSpaceDimension();
5405 int nElmtPntPrevious = 0;
5406 int nElmtCoefPrevious = 0;
5408 int nElmtPnt, nElmtCoef;
5409 for (
int nelmt = 0; nelmt < nTotElmt; nelmt++)
5411 nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
5412 nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
5415 if (nElmtPntPrevious != nElmtPnt || nElmtCoefPrevious != nElmtCoef ||
5416 (ElmtTypeNow != ElmtTypePrevious))
5418 if (nElmtPntPrevious != nElmtPnt)
5420 for (
int ndir = 0; ndir < nspacedim; ndir++)
5427 stdExp = (*m_exp)[nelmt]->GetStdExp();
5429 stdExp->DetShapeType(), *stdExp);
5431 ArrayStdMat[0] = stdExp->GetStdMatrix(matkey);
5432 ArrayStdMat_data[0] = ArrayStdMat[0]->GetPtr();
5439 ArrayStdMat[1] = stdExp->GetStdMatrix(matkey);
5440 ArrayStdMat_data[1] = ArrayStdMat[1]->GetPtr();
5445 stdExp->DetShapeType(),
5448 ArrayStdMat[2] = stdExp->GetStdMatrix(matkey);
5449 ArrayStdMat_data[2] = ArrayStdMat[2]->GetPtr();
5453 ElmtTypePrevious = ElmtTypeNow;
5454 nElmtPntPrevious = nElmtPnt;
5455 nElmtCoefPrevious = nElmtCoef;
5459 for (
int ndir = 0; ndir < nspacedim; ndir++)
5465 for (
int ndir = 0; ndir < nspacedim; ndir++)
5467 (*m_exp)[nelmt]->AlignVectorToCollapsedDir(
5468 ndir, inarray[ndir][nelmt], tmppnts);
5469 for (
int n = 0; n < nspacedim; n++)
5471 Vmath::Vadd(nElmtPnt, tmppnts[n], 1, projectedpnts[n], 1,
5472 projectedpnts[n], 1);
5476 for (
int ndir = 0; ndir < nspacedim; ndir++)
5479 (*m_exp)[nelmt]->MultiplyByQuadratureMetric(projectedpnts[ndir],
5480 projectedpnts[ndir]);
5483 for (
int np = 0; np < nElmtPnt; np++)
5485 NekDouble factor = projectedpnts[ndir][np];
5486 clmnArray = MatDataArray + np * nElmtCoef;
5487 clmnStdMatArray = ArrayStdMat_data[ndir] + np * nElmtCoef;
5488 Vmath::Svtvp(nElmtCoef, factor, clmnStdMatArray, 1, clmnArray,
5501 int nElmtPntPrevious = 0;
5502 int nElmtCoefPrevious = 0;
5503 int nTotElmt = (*m_exp).size();
5504 int nElmtPnt = (*m_exp)[0]->GetTotPoints();
5505 int nElmtCoef = (*m_exp)[0]->GetNcoeffs();
5511 for (
int nelmt = 0; nelmt < nTotElmt; nelmt++)
5513 nElmtCoef = (*m_exp)[nelmt]->GetNcoeffs();
5514 nElmtPnt = (*m_exp)[nelmt]->GetTotPoints();
5517 if (nElmtPntPrevious != nElmtPnt || nElmtCoefPrevious != nElmtCoef ||
5518 (ElmtTypeNow != ElmtTypePrevious))
5521 stdExp = (*m_exp)[nelmt]->GetStdExp();
5523 stdExp->DetShapeType(), *stdExp);
5526 stdMat_data = BwdMat->GetPtr();
5528 if (nElmtPntPrevious != nElmtPnt)
5533 ElmtTypePrevious = ElmtTypeNow;
5534 nElmtPntPrevious = nElmtPnt;
5535 nElmtCoefPrevious = nElmtCoef;
5538 (*m_exp)[nelmt]->MultiplyByQuadratureMetric(
5544 for (
int np = 0; np < nElmtPnt; np++)
5547 clmnArray = MatDataArray + np * nElmtCoef;
5548 clmnStdMatArray = stdMat_data + np * nElmtCoef;
5549 Vmath::Smul(nElmtCoef, factor, clmnStdMatArray, 1, clmnArray, 1);
5573 std::shared_ptr<LocalRegions::ExpansionVector> traceExp =
5574 tracelist->GetExp();
5575 int ntotTrace = (*traceExp).size();
5576 int nTracePnt, nTraceCoef;
5578 std::shared_ptr<LocalRegions::ExpansionVector> fieldExp =
GetExp();
5584 locTraceToTraceMap->GetLeftRightAdjacentExpId();
5586 locTraceToTraceMap->GetLeftRightAdjacentExpFlag();
5589 locTraceToTraceMap->GetTraceCoeffToLeftRightExpCoeffMap();
5591 locTraceToTraceMap->GetTraceCoeffToLeftRightExpCoeffSign();
5596 int MatIndex, nPnts;
5599 int nTracePntsTtl = tracelist->GetTotPoints();
5600 int nlocTracePts = locTraceToTraceMap->GetNLocTracePts();
5601 int nlocTracePtsFwd = locTraceToTraceMap->GetNFwdLocTracePts();
5602 int nlocTracePtsBwd = nlocTracePts - nlocTracePtsFwd;
5605 nlocTracePtsLR[0] = nlocTracePtsFwd;
5606 nlocTracePtsLR[1] = nlocTracePtsBwd;
5608 size_t nFwdBwdNonZero = 0;
5610 for (
int i = 0; i < 2; ++i)
5612 if (nlocTracePtsLR[i] > 0)
5614 tmpIndex[nFwdBwdNonZero] = i;
5620 for (
int i = 0; i < nFwdBwdNonZero; ++i)
5622 nlocTracePtsNonZeroIndex[i] = tmpIndex[i];
5628 for (
int k = 0; k < 2; ++k)
5633 for (
int k = 0; k < nFwdBwdNonZero; ++k)
5635 size_t i = nlocTracePtsNonZeroIndex[k];
5639 int nNumbElmt = fieldMat.size();
5642 for (
int i = 0; i < nNumbElmt; i++)
5644 ElmtMatDataArray[i] = fieldMat[i]->GetPtr();
5648 int nTraceCoefMax = 0;
5649 int nTraceCoefMin = std::numeric_limits<int>::max();
5655 for (
int nt = 0; nt < ntotTrace; nt++)
5657 nTraceCoef = (*traceExp)[nt]->GetNcoeffs();
5658 nTracePnt = tracelist->GetTotPoints(nt);
5659 int noffset = tracelist->GetPhys_Offset(nt);
5660 TraceCoefArray[nt] = nTraceCoef;
5661 TracePntArray[nt] = nTracePnt;
5662 TraceOffArray[nt] = noffset;
5663 FwdMatData[nt] = FwdMat[nt]->GetPtr();
5664 BwdMatData[nt] = BwdMat[nt]->GetPtr();
5665 if (nTraceCoef > nTraceCoefMax)
5667 nTraceCoefMax = nTraceCoef;
5669 if (nTraceCoef < nTraceCoefMin)
5671 nTraceCoefMin = nTraceCoef;
5674 WARNINGL1(nTraceCoefMax == nTraceCoefMin,
5675 "nTraceCoefMax!=nTraceCoefMin: Effeciency may be low ");
5677 int traceID, nfieldPnts, ElmtId, noffset;
5679 locTraceToTraceMap->GetLocTracephysToTraceIDMap();
5681 locTraceToTraceMap->GetLocTraceToFieldMap();
5684 for (
int k = 0; k < nFwdBwdNonZero; ++k)
5686 size_t i = nlocTracePtsNonZeroIndex[k];
5688 Vmath::Vcopy(nlocTracePtsLR[i], &fieldToLocTraceMap[0] + noffset, 1,
5689 &fieldToLocTraceMapLR[i][0], 1);
5690 noffset += nlocTracePtsLR[i];
5694 for (
int k = 0; k < nFwdBwdNonZero; ++k)
5696 size_t nlr = nlocTracePtsNonZeroIndex[k];
5698 for (
int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
5700 traceID = LocTracephysToTraceIDMap[nlr][nloc];
5701 nTraceCoef = TraceCoefArray[traceID];
5702 ElmtId = LRAdjExpid[nlr][traceID];
5704 nElmtCoef = ElmtCoefArray[ElmtId];
5705 nfieldPnts = fieldToLocTraceMapLR[nlr][nloc];
5706 nPnts = nfieldPnts - noffset;
5708 MatIndexArray[nlr][nloc] = nPnts * nElmtCoef;
5712 for (
int nc = 0; nc < nTraceCoefMin; nc++)
5714 for (
int nt = 0; nt < ntotTrace; nt++)
5716 nTraceCoef = TraceCoefArray[nt];
5717 nTracePnt = TracePntArray[nt];
5718 noffset = TraceOffArray[nt];
5719 Vmath::Vcopy(nTracePnt, &FwdMatData[nt][nc], nTraceCoef,
5720 &TraceFwdPhy[noffset], 1);
5721 Vmath::Vcopy(nTracePnt, &BwdMatData[nt][nc], nTraceCoef,
5722 &TraceBwdPhy[noffset], 1);
5725 for (
int k = 0; k < nFwdBwdNonZero; ++k)
5727 size_t i = nlocTracePtsNonZeroIndex[k];
5728 Vmath::Zero(nlocTracePtsLR[i], tmplocTrace[i], 1);
5734 for (
int k = 0; k < nFwdBwdNonZero; ++k)
5736 size_t nlr = nlocTracePtsNonZeroIndex[k];
5737 for (
int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
5739 traceID = LocTracephysToTraceIDMap[nlr][nloc];
5740 nTraceCoef = TraceCoefArray[traceID];
5741 ElmtId = LRAdjExpid[nlr][traceID];
5742 nrwAdjExp = elmtLRMap[nlr][traceID][nc];
5743 sign = elmtLRSign[nlr][traceID][nc];
5744 MatIndex = MatIndexArray[nlr][nloc] + nrwAdjExp;
5746 ElmtMatDataArray[ElmtId][MatIndex] -=
5747 sign * tmplocTrace[nlr][nloc];
5752 for (
int nc = nTraceCoefMin; nc < nTraceCoefMax; nc++)
5754 for (
int nt = 0; nt < ntotTrace; nt++)
5756 nTraceCoef = TraceCoefArray[nt];
5757 nTracePnt = TracePntArray[nt];
5758 noffset = TraceOffArray[nt];
5759 if (nc < nTraceCoef)
5761 Vmath::Vcopy(nTracePnt, &FwdMatData[nt][nc], nTraceCoef,
5762 &TraceFwdPhy[noffset], 1);
5763 Vmath::Vcopy(nTracePnt, &BwdMatData[nt][nc], nTraceCoef,
5764 &TraceBwdPhy[noffset], 1);
5773 for (
int k = 0; k < nFwdBwdNonZero; ++k)
5775 size_t i = nlocTracePtsNonZeroIndex[k];
5776 Vmath::Zero(nlocTracePtsLR[i], tmplocTrace[i], 1);
5781 for (
int k = 0; k < nFwdBwdNonZero; ++k)
5783 size_t nlr = nlocTracePtsNonZeroIndex[k];
5784 for (
int nloc = 0; nloc < nlocTracePtsLR[nlr]; nloc++)
5786 traceID = LocTracephysToTraceIDMap[nlr][nloc];
5787 nTraceCoef = TraceCoefArray[traceID];
5788 if (nc < nTraceCoef)
5790 ElmtId = LRAdjExpid[nlr][traceID];
5791 nrwAdjExp = elmtLRMap[nlr][traceID][nc];
5792 sign = -elmtLRSign[nlr][traceID][nc];
5793 MatIndex = MatIndexArray[nlr][nloc] + nrwAdjExp;
5795 ElmtMatDataArray[ElmtId][MatIndex] +=
5796 sign * tmplocTrace[nlr][nloc];
5808 int nelmtcoef, nelmtpnts, nelmtcoef0, nelmtpnts0;
5821 nelmtcoef0 = nelmtcoef;
5822 nelmtpnts0 = nelmtpnts;
5824 for (nelmt = 0; nelmt < (*m_exp).size(); ++nelmt)
5829 tmpMatQ = ElmtJacQuad[nelmt];
5830 tmpMatC = ElmtJacCoef[nelmt];
5832 MatQ_data = tmpMatQ->GetPtr();
5833 MatC_data = tmpMatC->GetPtr();
5835 if (nelmtcoef != nelmtcoef0)
5838 nelmtcoef0 = nelmtcoef;
5841 if (nelmtpnts != nelmtpnts0)
5844 nelmtpnts0 = nelmtpnts;
5847 for (
int np = 0; np < nelmtcoef; np++)
5849 Vmath::Vcopy(nelmtpnts, &MatQ_data[0] + np, nelmtcoef, &innarray[0],
5851 (*m_exp)[nelmt]->DivideByQuadratureMetric(innarray, innarray);
5852 (*m_exp)[nelmt]->IProductWRTDerivBase(dir, innarray, outarray);
5854 Vmath::Vadd(nelmtcoef, &outarray[0], 1, &MatC_data[0] + np,
5855 nelmtcoef, &MatC_data[0] + np, nelmtcoef);
5865 int nelmtcoef, nelmtpnts, nelmtcoef0, nelmtpnts0;
5878 nelmtcoef0 = nelmtcoef;
5879 nelmtpnts0 = nelmtpnts;
5881 for (nelmt = 0; nelmt < (*m_exp).size(); ++nelmt)
5886 tmpMatQ = ElmtJacQuad[nelmt];
5887 tmpMatC = ElmtJacCoef[nelmt];
5889 MatQ_data = tmpMatQ->GetPtr();
5890 MatC_data = tmpMatC->GetPtr();
5892 if (nelmtcoef != nelmtcoef0)
5895 nelmtcoef0 = nelmtcoef;
5898 if (nelmtpnts != nelmtpnts0)
5901 nelmtpnts0 = nelmtpnts;
5904 for (
int np = 0; np < nelmtcoef; np++)
5906 Vmath::Vcopy(nelmtpnts, &MatQ_data[0] + np, nelmtcoef, &innarray[0],
5908 (*m_exp)[nelmt]->DivideByQuadratureMetric(innarray, innarray);
5909 (*m_exp)[nelmt]->IProductWRTBase(innarray, outarray);
5911 Vmath::Vadd(nelmtcoef, &outarray[0], 1, &MatC_data[0] + np,
5912 nelmtcoef, &MatC_data[0] + np, nelmtcoef);
5932 int pt0 = (*m_exp)[i]->GetNumPoints(0);
5933 int pt1 = (*m_exp)[i]->GetNumPoints(1);
5934 int npt0 = (int)pt0 * scale;
5935 int npt1 = (int)pt1 * scale;
5938 npt0, (*
m_exp)[i]->GetPointsType(0));
5940 npt1, (*
m_exp)[i]->GetPointsType(1));
5944 newPointsKey0, newPointsKey1, &inarray[cnt],
5945 (*
m_exp)[i]->GetBasis(0)->GetPointsKey(),
5946 (*
m_exp)[i]->GetBasis(1)->GetPointsKey(), &outarray[cnt1]);
5958 int pt0 = (*m_exp)[i]->GetNumPoints(0);
5959 int pt1 = (*m_exp)[i]->GetNumPoints(1);
5960 int pt2 = (*m_exp)[i]->GetNumPoints(2);
5961 int npt0 = (int)pt0 * scale;
5962 int npt1 = (int)pt1 * scale;
5963 int npt2 = (int)pt2 * scale;
5966 npt0, (*
m_exp)[i]->GetPointsType(0));
5968 npt1, (*
m_exp)[i]->GetPointsType(1));
5970 npt2, (*
m_exp)[i]->GetPointsType(2));
5974 newPointsKey0, newPointsKey1, newPointsKey2, &inarray[cnt],
5975 (*
m_exp)[i]->GetBasis(0)->GetPointsKey(),
5976 (*
m_exp)[i]->GetBasis(1)->GetPointsKey(),
5977 (*
m_exp)[i]->GetBasis(2)->GetPointsKey(), &outarray[cnt1]);
5979 cnt += npt0 * npt1 * npt2;
5980 cnt1 += pt0 * pt1 * pt2;
#define WARNINGL1(condition, msg)
#define ASSERTL0(condition, msg)
#define NEKERROR(type, msg)
Assert Level 0 – Fundamental assert which is used whether in FULLDEBUG, DEBUG or OPT compilation mode...
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
#define LIKWID_MARKER_START(regionTag)
#define LIKWID_MARKER_STOP(regionTag)
#define sign(a, b)
return the sign(b)*a
COLLECTIONS_EXPORT OperatorImpMap SetWithTimings(std::vector< StdRegions::StdExpansionSharedPtr > pGeom, OperatorImpMap &impTypes, bool verbose=true)
unsigned int GetMaxCollectionSize()
COLLECTIONS_EXPORT void UpdateOptFile(std::string sessName, LibUtilities::CommSharedPtr &comm)
COLLECTIONS_EXPORT OperatorImpMap GetOperatorImpMap(StdRegions::StdExpansionSharedPtr pExp)
Get Operator Implementation Map from XMl or using default;.
Describes the specification for a Basis.
int GetNumPoints() const
Return points order at which basis is defined.
BasisType GetBasisType() const
Return type of expansion basis.
PointsKey GetPointsKey() const
Return distribution of points.
int GetNumModes() const
Returns the order of the basis.
static const std::string GetFileType(const std::string &filename, CommSharedPtr comm)
Determine file type of given input file.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
Defines a specification for a set of points.
void AccumulateRegion(std::string, int iolevel=0)
Accumulate elapsed time for a region.
General purpose memory allocation routines with the ability to allocate from thread specific memory p...
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
Base class for all multi-elemental spectral/hp expansions.
virtual void v_Reset()
Reset geometry information, metrics, matrix managers and geometry information.
Array< OneD, NekDouble > m_coeffs
Concatenation of all local expansion coefficients.
void AddTraceJacToElmtJac(const Array< OneD, const DNekMatSharedPtr > &FwdMat, const Array< OneD, const DNekMatSharedPtr > &BwdMat, Array< OneD, DNekMatSharedPtr > &fieldMat)
inverse process of v_GetFwdBwdTracePhys. Given Trace integration of Fwd and Bwd Jacobian,...
void IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to all local expansion modes .
const Array< OneD, const std::shared_ptr< ExpList > > & GetBndCondExpansions()
std::shared_ptr< ExpList > & GetTrace()
int GetNcoeffs(void) const
Returns the total number of local degrees of freedom .
virtual void v_ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
virtual 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)
const DNekScalBlkMatSharedPtr & GetBlockMatrix(const GlobalMatrixKey &gkey)
virtual void v_AddFwdBwdTraceIntegral(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &outarray)
ExpList(const ExpansionType Type=eNoType)
The default constructor using a type.
void ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
Extract the data in fielddata into the coeffs.
void CreateCollections(Collections::ImplementationType ImpType=Collections::eNoImpType)
Construct collections of elements containing a single element type and polynomial order from the list...
virtual void v_GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
void ExponentialFilter(Array< OneD, NekDouble > &array, const NekDouble alpha, const NekDouble exponent, const NekDouble cutoff)
virtual void v_HomogeneousBwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
std::shared_ptr< GlobalMatrix > GenGlobalMatrix(const GlobalMatrixKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
Generates a global matrix from the given key and map.
virtual void v_HomogeneousFwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, bool Shuff=true, bool UnShuff=true)
virtual void v_IProductWRTBase(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void AddRightIPTPhysDerivBase(const int dir, const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
virtual void v_ExtractTracePhys(Array< OneD, NekDouble > &outarray)
int GetExpIndex(const Array< OneD, const NekDouble > &gloCoord, NekDouble tol=0.0, bool returnNearestElmt=false, int cachedId=-1, NekDouble maxDistance=1e6)
This function returns the index of the local elemental expansion containing the arbitrary point given...
virtual void v_WriteTecplotConnectivity(std::ostream &outfile, int expansion)
BlockMatrixMapShPtr m_blockMat
virtual Array< OneD, const unsigned int > v_GetYIDs(void)
void MultiplyByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
multiply the metric jacobi and quadrature weights
void GetBoundaryToElmtMap(Array< OneD, int > &ElmtID, Array< OneD, int > &EdgeID)
void WriteVtkPieceFooter(std::ostream &outfile, int expansion)
virtual void v_ExtractPhysToBnd(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bnd)
void GetBwdWeight(Array< OneD, NekDouble > &weightAver, Array< OneD, NekDouble > &weightJump)
Get the weight value for boundary conditions for boundary average and jump calculations.
virtual void v_AddTraceIntegralToOffDiag(const Array< OneD, const NekDouble > &FwdFlux, const Array< OneD, const NekDouble > &BwdFlux, Array< OneD, NekDouble > &outarray)
virtual std::vector< bool > & v_GetLeftAdjacentTraces(void)
void InitialiseExpVector(const SpatialDomains::ExpansionInfoMap &expmap)
Define a list of elements using the geometry and basis key information in expmap;.
void SetupCoeffPhys(bool DeclareCoeffPhysArrays=true, bool SetupOffsets=true)
Definition of the total number of degrees of freedom and quadrature points and offsets to access data...
static SpatialDomains::BoundaryConditionShPtr GetBoundaryCondition(const SpatialDomains::BoundaryConditionCollection &collection, unsigned int index, const std::string &variable)
virtual void v_DealiasedDotProd(const Array< OneD, Array< OneD, NekDouble >> &inarray1, const Array< OneD, Array< OneD, NekDouble >> &inarray2, Array< OneD, Array< OneD, NekDouble >> &outarray)
virtual void v_CurlCurl(Array< OneD, Array< OneD, NekDouble >> &Vel, Array< OneD, Array< OneD, NekDouble >> &Q)
virtual void v_PhysDirectionalDeriv(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_GetBCValues(Array< OneD, NekDouble > &BndVals, const Array< OneD, NekDouble > &TotField, int BndID)
virtual void v_GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
virtual void v_WriteTecplotHeader(std::ostream &outfile, std::string var="")
virtual std::shared_ptr< ExpList > & v_GetPlane(int n)
int GetCoeff_Offset(int n) const
Get the start offset position for a global list of m_coeffs correspoinding to element n.
int GetExpSize(void)
This function returns the number of elements in the expansion.
virtual void v_AddTraceQuadPhysToOffDiag(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
virtual void v_ClearGlobalLinSysManager(void)
void GeneralGetFieldDefinitions(std::vector< LibUtilities::FieldDefinitionsSharedPtr > &fielddef, int NumHomoDir=0, Array< OneD, LibUtilities::BasisSharedPtr > &HomoBasis=LibUtilities::NullBasisSharedPtr1DArray, std::vector< NekDouble > &HomoLen=LibUtilities::NullNekDoubleVector, bool homoStrips=false, std::vector< unsigned int > &HomoSIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoZIDs=LibUtilities::NullUnsignedIntVector, std::vector< unsigned int > &HomoYIDs=LibUtilities::NullUnsignedIntVector)
virtual void v_GetPeriodicEntities(PeriodicMap &periodicVerts, PeriodicMap &periodicEdges, PeriodicMap &periodicFaces)
void IProductWRTDerivBase(const int dir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the inner product of a function with respect to the derivative (in directio...
virtual void v_LocalToGlobal(bool UseComm)
int GetNumElmts(void)
This function returns the number of elements in the expansion which may be different for a homogeoeno...
virtual const Array< OneD, const SpatialDomains::BoundaryConditionShPtr > & v_GetBndConditions()
virtual ~ExpList()
The default destructor.
virtual void v_FillBndCondFromField()
virtual void v_FwdTransLocalElmt(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_ImposeDirichletConditions(Array< OneD, NekDouble > &outarray)
virtual void v_EvaluateBoundaryConditions(const NekDouble time=0.0, const std::string varName="", const NekDouble x2_in=NekConstants::kNekUnsetDouble, const NekDouble x3_in=NekConstants::kNekUnsetDouble)
void ClearGlobalLinSysManager(void)
virtual const std::vector< bool > & v_GetLeftAdjacentFaces(void) const
virtual void v_GetNormals(Array< OneD, Array< OneD, NekDouble >> &normals)
Populate normals with the normals of all expansions.
Array< OneD, int > m_coeff_offset
Offset of elemental data into the array m_coeffs.
virtual NekDouble v_VectorFlux(const Array< OneD, Array< OneD, NekDouble >> &inarray)
std::shared_ptr< AssemblyMapDG > & GetTraceMap(void)
virtual std::shared_ptr< AssemblyMapDG > & v_GetTraceMap()
Array< OneD, std::pair< int, int > > m_coeffsToElmt
m_coeffs to elemental value map
virtual void v_GetFwdBwdTracePhys(Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
NekDouble PhysEvaluate(const Array< OneD, const NekDouble > &coords, const Array< OneD, const NekDouble > &phys)
virtual NekDouble v_Integral(const Array< OneD, const NekDouble > &inarray)
virtual std::map< int, RobinBCInfoSharedPtr > v_GetRobinBCInfo(void)
virtual const Array< OneD, const int > & v_GetTraceBndMap()
virtual Array< OneD, const NekDouble > v_HomogeneousEnergy(void)
virtual NekDouble v_GetHomoLen(void)
void ExtractFileBCs(const std::string &fileName, LibUtilities::CommSharedPtr comm, const std::string &varName, const std::shared_ptr< ExpList > locExpList)
const LocTraceToTraceMapSharedPtr & GetLocTraceToTraceMap() const
virtual void v_LinearAdvectionReactionSolve(const Array< OneD, Array< OneD, NekDouble >> &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
virtual Array< OneD, const unsigned int > v_GetZIDs(void)
void ExtractCoeffsToCoeffs(const std::shared_ptr< ExpList > &fromExpList, const Array< OneD, const NekDouble > &fromCoeffs, Array< OneD, NekDouble > &toCoeffs)
Extract the data from fromField using fromExpList the coeffs using the basic ExpList Elemental expans...
void GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1=NullNekDouble1DArray, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
This function calculates the coordinates of all the elemental quadrature points .
std::shared_ptr< GlobalLinSys > GenGlobalBndLinSys(const GlobalLinSysKey &mkey, const AssemblyMapSharedPtr &locToGloMap)
Generate a GlobalLinSys from information provided by the key "mkey" and the mapping provided in LocTo...
virtual void v_SetBndCondBwdWeight(const int index, const NekDouble value)
void GetMatIpwrtDeriveBase(const Array< OneD, const Array< OneD, NekDouble >> &inarray, const int nDirctn, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
virtual void v_SmoothField(Array< OneD, NekDouble > &field)
virtual void v_AppendFieldData(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata)
std::vector< bool > m_collectionsDoInit
Vector of bools to act as an initialise on first call flag.
bool m_physState
The state of the array m_phys.
void AddRightIPTBaseMatrix(const Array< OneD, const DNekMatSharedPtr > ElmtJacQuad, Array< OneD, DNekMatSharedPtr > ElmtJacCoef)
void GetLocTraceFromTracePts(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &locTraceFwd, Array< OneD, NekDouble > &locTraceBwd)
virtual void v_WriteVtkPieceData(std::ostream &outfile, int expansion, std::string var)
virtual void v_WriteVtkPieceHeader(std::ostream &outfile, int expansion, int istrip)
int GetShapeDimension()
This function returns the dimension of the shape of the element eid.
NekDouble H1(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
Calculates the error of the global spectral/hp element approximation.
std::vector< int > m_coll_phys_offset
Offset of elemental data into the array m_phys.
LibUtilities::CommSharedPtr m_comm
Communicator.
std::shared_ptr< GlobalLinSys > GenGlobalLinSys(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
This operation constructs the global linear system of type mkey.
void WriteVtkHeader(std::ostream &outfile)
std::shared_ptr< LocalRegions::ExpansionVector > m_exp
The list of local expansions.
virtual void v_GetCoords(Array< OneD, NekDouble > &coord_0, Array< OneD, NekDouble > &coord_1, Array< OneD, NekDouble > &coord_2=NullNekDouble1DArray)
int m_ncoeffs
The total number of local degrees of freedom. m_ncoeffs .
void GenerateElementVector(const int ElementID, const NekDouble scalar1, const NekDouble scalar2, Array< OneD, NekDouble > &outarray)
Generate vector v such that v[i] = scalar1 if i is in the element < ElementID. Otherwise,...
std::shared_ptr< DNekMat > GenGlobalMatrixFull(const GlobalLinSysKey &mkey, const std::shared_ptr< AssemblyMapCG > &locToGloMap)
virtual void v_NormVectorIProductWRTBase(Array< OneD, const NekDouble > &V1, Array< OneD, const NekDouble > &V2, Array< OneD, NekDouble > &outarray, int BndID)
virtual void v_ExtractElmtToBndPhys(const int i, const Array< OneD, NekDouble > &elmt, Array< OneD, NekDouble > &boundary)
virtual void v_ExtractPhysToBndElmt(const int i, const Array< OneD, const NekDouble > &phys, Array< OneD, NekDouble > &bndElmt)
virtual void v_AddTraceIntegral(const Array< OneD, const NekDouble > &Fx, const Array< OneD, const NekDouble > &Fy, Array< OneD, NekDouble > &outarray)
virtual void v_FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
virtual void v_PhysGalerkinProjection1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
void GeneralMatrixOp(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function calculates the result of the multiplication of a matrix of type specified by mkey with ...
virtual void v_SetUpPhysNormals()
: Set up a normal along the trace elements between two elements at elemental level
virtual void v_GlobalToLocal(void)
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void)
const std::shared_ptr< LocalRegions::ExpansionVector > GetExp() const
This function returns the vector of elements in the expansion.
virtual void v_HelmSolve(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const StdRegions::ConstFactorMap &factors, const StdRegions::VarCoeffMap &varcoeff, const MultiRegions::VarFactorsMap &varfactors, const Array< OneD, const NekDouble > &dirForcing, const bool PhysSpaceForcing)
SpatialDomains::MeshGraphSharedPtr m_graph
Mesh associated with this expansion list.
virtual void v_MultiplyByInvMassMatrix(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_AddTraceQuadPhysToField(const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &field)
virtual const Array< OneD, const NekDouble > & v_GetBndCondBwdWeight()
virtual void v_GetMovingFrames(const SpatialDomains::GeomMMF MMFdir, const Array< OneD, const NekDouble > &CircCentre, Array< OneD, Array< OneD, NekDouble >> &outarray)
virtual NekDouble v_L2(const Array< OneD, const NekDouble > &phys, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
void WriteVtkFooter(std::ostream &outfile)
void MultiplyByBlockMatrix(const GlobalMatrixKey &gkey, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual std::shared_ptr< ExpList > & v_GetTrace()
void GetElmtNormalLength(Array< OneD, NekDouble > &lengthsFwd, Array< OneD, NekDouble > &lengthsBwd)
Get the length of elements in boundary normal direction.
virtual std::vector< LibUtilities::FieldDefinitionsSharedPtr > v_GetFieldDefinitions(void)
std::vector< int > m_coll_coeff_offset
Offset of elemental data into the array m_coeffs.
void IProductWRTDirectionalDerivBase(const Array< OneD, const NekDouble > &direction, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Directional derivative along a given direction.
void FillBwdWithBwdWeight(Array< OneD, NekDouble > &weightave, Array< OneD, NekDouble > &weightjmp)
Fill Bwd with boundary conditions.
virtual void v_BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual Array< OneD, SpatialDomains::BoundaryConditionShPtr > & v_UpdateBndConditions()
virtual void v_LinearAdvectionDiffusionReactionSolve(const Array< OneD, Array< OneD, NekDouble >> &velocity, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray, const NekDouble lambda, const Array< OneD, const NekDouble > &dirForcing=NullNekDouble1DArray)
virtual void v_WriteTecplotField(std::ostream &outfile, int expansion)
std::shared_ptr< ExpList > GetSharedThisPtr()
Returns a shared pointer to the current object.
Collections::CollectionVector m_collections
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array< OneD, NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual void v_GetBoundaryNormals(int i, Array< OneD, Array< OneD, NekDouble >> &normals)
void ApplyGeomInfo()
Apply geometry information to each expansion.
Array< OneD, const unsigned int > GetZIDs(void)
This function returns a vector containing the wave numbers in z-direction associated with the 3D homo...
virtual void v_PhysDeriv(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d0, Array< OneD, NekDouble > &out_d1, Array< OneD, NekDouble > &out_d2)
virtual void v_FwdTransBndConstrained(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
std::unordered_map< int, int > m_elmtToExpId
Mapping from geometry ID of element to index inside m_exp.
void GetDiagMatIpwrtBase(const Array< OneD, const Array< OneD, NekDouble >> &inarray, Array< OneD, DNekMatSharedPtr > &mtxPerVar)
virtual void v_WriteTecplotZone(std::ostream &outfile, int expansion)
LibUtilities::SessionReaderSharedPtr m_session
Session.
virtual void v_PeriodicBwdCopy(const Array< OneD, const NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd)
std::map< int, RobinBCInfoSharedPtr > GetRobinBCInfo()
virtual void v_FillBwdWithBoundCond(const Array< OneD, NekDouble > &Fwd, Array< OneD, NekDouble > &Bwd, bool PutFwdInBwdOnBCs)
void DivideByQuadratureMetric(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
Divided by the metric jacobi and quadrature weights.
virtual const Array< OneD, const std::shared_ptr< ExpList > > & v_GetBndCondExpansions(void)
Array< OneD, int > m_phys_offset
Offset of elemental data into the array m_phys.
virtual void v_SetHomoLen(const NekDouble lhom)
void Upwind(const Array< OneD, const NekDouble > &Vn, const Array< OneD, const NekDouble > &Fwd, const Array< OneD, const NekDouble > &Bwd, Array< OneD, NekDouble > &Upwind)
const DNekScalBlkMatSharedPtr GenBlockMatrix(const GlobalMatrixKey &gkey)
This function assembles the block diagonal matrix of local matrices of the type mtype.
void PhysDeriv(Direction edir, const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &out_d)
void MultiplyByElmtInvMass(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally mulplies the coefficient space of Sin my the elemental inverse of the mass ...
void BwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
This function elementally evaluates the backward transformation of the global spectral/hp element exp...
virtual void v_GetBndElmtExpansion(int i, std::shared_ptr< ExpList > &result, const bool DeclareCoeffPhysArrays)
virtual void v_FwdTrans(const Array< OneD, const NekDouble > &inarray, Array< OneD, NekDouble > &outarray)
virtual std::shared_ptr< ExpList > & v_UpdateBndCondExpansion(int i)
int GetPhys_Offset(int n) const
Get the start offset position for a global list of m_phys correspoinding to element n.
virtual const std::shared_ptr< LocTraceToTraceMap > & v_GetLocTraceToTraceMap(void) const
virtual void v_ExtractDataToCoeffs(LibUtilities::FieldDefinitionsSharedPtr &fielddef, std::vector< NekDouble > &fielddata, std::string &field, Array< OneD, NekDouble > &coeffs)
Extract data from raw field data into expansion list.
ExpansionType m_expType
Exapnsion type.
Array< OneD, NekDouble > m_phys
The global expansion evaluated at the quadrature points.
ExpansionType GetExpType(void)
Returns the type of the expansion.
NekDouble Linf(const Array< OneD, const NekDouble > &inarray, const Array< OneD, const NekDouble > &soln=NullNekDouble1DArray)
This function calculates the error of the global spectral/hp element approximation.
virtual void v_DealiasedProd(const Array< OneD, NekDouble > &inarray1, const Array< OneD, NekDouble > &inarray2, Array< OneD, NekDouble > &outarray)
int GetTotPoints(void) const
Returns the total number of quadrature points m_npoints .
int GetCoordim(int eid)
This function returns the dimension of the coordinates of the element eid.
Describe a linear system.
GlobalSysSolnType GetGlobalSysSolnType() const
Return the associated solution type.
Describes a matrix with ordering defined by a local to global map.
const StdRegions::ConstFactorMap & GetConstFactors() const
Returns all the constants.
int GetNVarCoeffs() const
LibUtilities::ShapeType GetShapeType() const
Return the expansion type associated with key.
const StdRegions::VarCoeffMap & GetVarCoeffs() const
StdRegions::MatrixType GetMatrixType() const
Return the matrix type.
static void Dscal(const int &n, const double &alpha, double *x, const int &incx)
BLAS level 1: x = alpha x.
std::map< OperatorType, ImplementationType > OperatorImpMap
void PhysGalerkinProject3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
std::shared_ptr< FieldIO > FieldIOSharedPtr
static const BasisKey NullBasisKey(eNoBasisType, 0, NullPointsKey)
Defines a null basis with no type or points.
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 ...
std::shared_ptr< SessionReader > SessionReaderSharedPtr
void Interp3D(const BasisKey &fbasis0, const BasisKey &fbasis1, const BasisKey &fbasis2, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, const BasisKey &tbasis2, Array< OneD, NekDouble > &to)
this function interpolates a 3D function evaluated at the quadrature points of the 3D basis,...
std::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
void Interp2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
this function interpolates a 2D function evaluated at the quadrature points of the 2D basis,...
int GetNumberOfCoefficients(ShapeType shape, std::vector< unsigned int > &modes, int offset=0)
FieldIOFactory & GetFieldIOFactory()
Returns the FieldIO factory.
std::shared_ptr< Transposition > TranspositionSharedPtr
@ eNodalTriElec
2D Nodal Electrostatic Points on a Triangle
std::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
void PhysGalerkinProject2D(const BasisKey &fbasis0, const BasisKey &fbasis1, const Array< OneD, const NekDouble > &from, const BasisKey &tbasis0, const BasisKey &tbasis1, Array< OneD, NekDouble > &to)
@ eGauss_Lagrange
Lagrange Polynomials using the Gauss points.
@ eOrtho_A
Principle Orthogonal Functions .
@ eGLL_Lagrange
Lagrange for SEM basis .
std::shared_ptr< PointExp > PointExpSharedPtr
std::shared_ptr< Expansion > ExpansionSharedPtr
std::shared_ptr< Expansion2D > Expansion2DSharedPtr
std::shared_ptr< Expansion0D > Expansion0DSharedPtr
std::shared_ptr< Expansion1D > Expansion1DSharedPtr
std::shared_ptr< Expansion3D > Expansion3DSharedPtr
std::vector< ExpansionSharedPtr > ExpansionVector
MultiRegions::Direction const DirCartesianMap[]
@ eSIZE_GlobalSysSolnType
void AlignFace(const StdRegions::Orientation orient, const int nquad1, const int nquad2, const Array< OneD, const NekDouble > &in, Array< OneD, NekDouble > &out)
Helper function to re-align face to a given orientation.
static ExpListSharedPtr NullExpListSharedPtr
std::shared_ptr< AssemblyMapCG > AssemblyMapCGSharedPtr
std::shared_ptr< RobinBCInfo > RobinBCInfoSharedPtr
const char *const GlobalSysSolnTypeMap[]
std::shared_ptr< GlobalLinSys > GlobalLinSysSharedPtr
Pointer to a GlobalLinSys object.
std::shared_ptr< GlobalMatrix > GlobalMatrixSharedPtr
Shared pointer to a GlobalMatrix object.
GlobalLinSysFactory & GetGlobalLinSysFactory()
static LocTraceToTraceMapSharedPtr NullLocTraceToTraceMapSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
std::shared_ptr< LocTraceToTraceMap > LocTraceToTraceMapSharedPtr
std::map< StdRegions::ConstFactorType, Array< OneD, NekDouble > > VarFactorsMap
std::shared_ptr< AssemblyMap > AssemblyMapSharedPtr
std::map< int, std::vector< PeriodicEntity > > PeriodicMap
std::map< GlobalMatrixKey, DNekScalBlkMatSharedPtr > BlockMatrixMap
A map between global matrix keys and their associated block matrices.
static const NekDouble kNekZeroTol
std::shared_ptr< QuadGeom > QuadGeomSharedPtr
GeomMMF
Principle direction for MMF.
std::shared_ptr< BoundaryConditionBase > BoundaryConditionShPtr
std::shared_ptr< PrismGeom > PrismGeomSharedPtr
std::map< int, CompositeSharedPtr > CompositeMap
std::map< int, BoundaryConditionMapShPtr > BoundaryConditionCollection
@ eDeformed
Geometry is curved or has non-constant factors.
std::shared_ptr< SegGeom > SegGeomSharedPtr
std::shared_ptr< PyrGeom > PyrGeomSharedPtr
std::shared_ptr< ExpansionInfo > ExpansionInfoShPtr
std::shared_ptr< TetGeom > TetGeomSharedPtr
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
std::shared_ptr< BoundaryConditionMap > BoundaryConditionMapShPtr
std::shared_ptr< HexGeom > HexGeomSharedPtr
std::shared_ptr< PointGeom > PointGeomSharedPtr
std::shared_ptr< Geometry2D > Geometry2DSharedPtr
std::shared_ptr< TriGeom > TriGeomSharedPtr
std::shared_ptr< Geometry1D > Geometry1DSharedPtr
std::map< int, ExpansionInfoShPtr > ExpansionInfoMap
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
std::map< StdRegions::VarCoeffType, Array< OneD, NekDouble > > VarCoeffMap
std::map< ConstFactorType, NekDouble > ConstFactorMap
@ eDir1BwdDir2_Dir2BwdDir1
@ eDir1BwdDir1_Dir2BwdDir2
@ eDir1BwdDir2_Dir2FwdDir1
@ eDir1FwdDir1_Dir2BwdDir2
@ eDir1BwdDir1_Dir2FwdDir2
@ eDir1FwdDir2_Dir2FwdDir1
@ eDir1FwdDir2_Dir2BwdDir1
The above copyright notice and this permission notice shall be included.
NekMatrix< NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag >, BlockMatrixTag > DNekScalBlkMat
NekMatrix< NekMatrix< NekDouble, StandardMatrixTag >, ScaledMatrixTag > DNekScalMat
std::shared_ptr< DNekScalMat > DNekScalMatSharedPtr
std::pair< IndexType, IndexType > CoordType
std::shared_ptr< DNekScalBlkMat > DNekScalBlkMatSharedPtr
@ ePOSITIVE_DEFINITE_SYMMETRIC_BANDED
@ ePOSITIVE_DEFINITE_SYMMETRIC
std::map< CoordType, NekDouble > COOMatType
static Array< OneD, NekDouble > NullNekDouble1DArray
std::shared_ptr< DNekMat > DNekMatSharedPtr
void Svtvp(int n, const T alpha, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
svtvp (scalar times vector plus vector): z = alpha*x + y
void Neg(int n, T *x, const int incx)
Negate x = -x.
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
void Scatr(int n, const T *x, const int *y, T *z)
Scatter vector z[y[i]] = x[i].
void Vadd(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Add vector z = x+y.
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
void Zero(int n, T *x, const int incx)
Zero vector.
T Vmax(int n, const T *x, const int incx)
Return the maximum element in x – called vmax to avoid conflict with max.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
scalarT< T > abs(scalarT< T > in)
scalarT< T > sqrt(scalarT< T > in)
Used to lookup the create function in NekManager.
std::shared_ptr< RobinBCInfo > next