Nektar++
DriverAdaptive.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: DriverAdaptive.cpp
4//
5// For more information, please see: http://www.nektar.info
6//
7// The MIT License
8//
9// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10// Department of Aeronautics, Imperial College London (UK), and Scientific
11// Computing and Imaging Institute, University of Utah (USA).
12//
13// Permission is hereby granted, free of charge, to any person obtaining a
14// copy of this software and associated documentation files (the "Software"),
15// to deal in the Software without restriction, including without limitation
16// the rights to use, copy, modify, merge, publish, distribute, sublicense,
17// and/or sell copies of the Software, and to permit persons to whom the
18// Software is furnished to do so, subject to the following conditions:
19//
20// The above copyright notice and this permission notice shall be included
21// in all copies or substantial portions of the Software.
22//
23// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29// DEALINGS IN THE SOFTWARE.
30//
31// Description: Driver class for adaptive solver
32//
33///////////////////////////////////////////////////////////////////////////////
34
35#include <iomanip>
36
47
48namespace Nektar::SolverUtils
49{
50
51std::string DriverAdaptive::className =
56
57/**
58 *
59 */
63 : Driver(pSession, pGraph)
64{
65}
66
67/**
68 *
69 */
70void DriverAdaptive::v_InitObject(std::ostream &out)
71{
73}
74
75/**
76 *
77 */
78void DriverAdaptive::v_Execute(std::ostream &out)
79{
81 NekDouble CPUtime;
82
83 m_equ[0]->PrintSummary(out);
84
85 // First run using original order
86 timer.Start();
87 m_equ[0]->DoInitialise();
88
89 // Obtain initial time in case a restart was used
90 NekDouble startTime = m_equ[0]->GetTime();
91 m_equ[0]->DoSolve();
92
93 // Load session parameters and solver info
94 bool isHomogeneous1D;
95 int nRuns, minP, maxP, sensorVar;
96 NekDouble lowerTol, upperTol;
97 m_session->LoadParameter("NumRuns", nRuns, 1);
98 m_session->LoadParameter("AdaptiveMinModes", minP, 4);
99 m_session->LoadParameter("AdaptiveMaxModes", maxP, 12);
100 m_session->LoadParameter("AdaptiveLowerTolerance", lowerTol, 1e-8);
101 m_session->LoadParameter("AdaptiveUpperTolerance", upperTol, 1e-6);
102 m_session->LoadParameter("AdaptiveSensorVariable", sensorVar, 0);
103 m_session->MatchSolverInfo("Homogeneous", "1D", isHomogeneous1D, false);
104
105 // Get number of elements and planes
106 int nExp, nPlanes;
107 if (isHomogeneous1D)
108 {
109 nExp = m_equ[0]->UpdateFields()[0]->GetPlane(0)->GetExpSize();
110 nPlanes = m_equ[0]->UpdateFields()[0]->GetZIDs().size();
111 }
112 else
113 {
114 nExp = m_equ[0]->UpdateFields()[0]->GetExpSize();
115 nPlanes = 1;
116 }
117 int expdim = m_equ[0]->UpdateFields()[0]->GetGraph()->GetMeshDimension();
118
119 int nFields = m_equ[0]->UpdateFields().size();
120 int numSteps = m_session->GetParameter("NumSteps");
121 NekDouble period = m_session->GetParameter("TimeStep") * numSteps;
122
123 Array<OneD, NekDouble> coeffs, phys, physReduced, tmpArray;
124
125 // Get mapping
127 mapping = GlobalMapping::Mapping::Load(m_session, m_equ[0]->UpdateFields());
128
129 // Adaptive loop
130 Array<OneD, int> P(expdim);
131 Array<OneD, int> numPoints(expdim);
134 for (int i = 1; i < nRuns; i++)
135 {
136 // Get field expansions
138 m_equ[0]->UpdateFields();
139
140 // Determine the change to be applied in the order
141 std::map<int, int> deltaP;
142 int offset = 0;
143 for (int n = 0; n < nExp; n++)
144 {
145 offset = fields[sensorVar]->GetPhys_Offset(n);
146 Exp = fields[sensorVar]->GetExp(n);
147
148 for (int k = 0; k < expdim; ++k)
149 {
150 P[k] = Exp->GetBasis(k)->GetNumModes();
151 numPoints[k] = Exp->GetBasis(k)->GetNumPoints();
152 ptsKey[k] = LibUtilities::PointsKey(
153 numPoints[k], Exp->GetBasis(k)->GetPointsType());
154 }
155
156 // Declare orthogonal basis.
158 switch (Exp->GetGeom()->GetShapeType())
159 {
161 {
163 ptsKey[0]);
165 ptsKey[1]);
166 OrthoExp = MemoryManager<
167 StdRegions::StdQuadExp>::AllocateSharedPtr(Ba, Bb);
168 break;
169 }
171 {
173 ptsKey[0]);
175 ptsKey[1]);
176 OrthoExp =
178 Ba, Bb);
179 break;
180 }
182 {
184 ptsKey[0]);
186 ptsKey[1]);
188 ptsKey[2]);
189 OrthoExp =
191 Ba, Bb, Bc);
192 break;
193 }
195 {
197 ptsKey[0]);
199 ptsKey[1]);
201 ptsKey[2]);
202 OrthoExp =
204 Ba, Bb, Bc);
205 break;
206 }
208 {
210 ptsKey[0]);
212 ptsKey[1]);
214 ptsKey[2]);
215 OrthoExp = MemoryManager<
216 StdRegions::StdPrismExp>::AllocateSharedPtr(Ba, Bb, Bc);
217 break;
218 }
220 {
222 ptsKey[0]);
224 ptsKey[1]);
226 ptsKey[2]);
227 OrthoExp =
229 Ba, Bb, Bc);
230 break;
231 }
232 default:
233 ASSERTL0(false, "Shape not supported.");
234 break;
235 }
236
237 int nq = OrthoExp->GetTotPoints();
238
239 NekDouble error = 0;
240 NekDouble fac = 0;
241 NekDouble tmp = 0;
242
243 coeffs = Array<OneD, NekDouble>(OrthoExp->GetNcoeffs());
244 physReduced = Array<OneD, NekDouble>(OrthoExp->GetTotPoints());
245 tmpArray = Array<OneD, NekDouble>(OrthoExp->GetTotPoints(), 0.0);
246
247 // Refinement based only on one variable
248 for (int plane = 0; plane < nPlanes; plane++)
249 {
250 if (isHomogeneous1D)
251 {
252 phys =
253 fields[sensorVar]->GetPlane(plane)->GetPhys() + offset;
254 }
255 else
256 {
257 phys = fields[sensorVar]->GetPhys() + offset;
258 }
259
260 // Project solution to lower order
261 OrthoExp->FwdTrans(phys, coeffs);
262 OrthoExp->BwdTrans(coeffs, physReduced);
263
264 // Calculate error =||phys-physReduced||^2 / ||phys||^2
265 Vmath::Vsub(nq, phys, 1, physReduced, 1, tmpArray, 1);
266 Vmath::Vmul(nq, tmpArray, 1, tmpArray, 1, tmpArray, 1);
267 tmp = Exp->Integral(tmpArray);
268
269 Vmath::Vmul(nq, phys, 1, phys, 1, tmpArray, 1);
270 fac = Exp->Integral(tmpArray);
271
272 tmp = abs(tmp / fac);
273
274 if (tmp != tmp)
275 {
276 ASSERTL0(false,
277 "Adaptive procedure encountered NaN value.");
278 }
279
280 // Get maximum value along planes
281 error = (tmp > error) ? tmp : error;
282 }
283
284 // Combine planes from different processes
285 m_session->GetComm()->GetColumnComm()->AllReduce(
287
288 // Override tolerances if function is defined
289 if (m_session->DefinesFunction("AdaptiveLowerTolerance"))
290 {
291 int nq = Exp->GetTotPoints();
292
293 // Obtain points from the the element
294 Array<OneD, NekDouble> xc0, xc1, xc2;
295 xc0 = Array<OneD, NekDouble>(nq, 0.0);
296 xc1 = Array<OneD, NekDouble>(nq, 0.0);
297 xc2 = Array<OneD, NekDouble>(nq, 0.0);
298 Exp->GetCoords(xc0, xc1, xc2);
299
300 // Evaluate function from session file
301 Array<OneD, NekDouble> tolerance(nq, 0.0);
303 m_session->GetFunction("AdaptiveLowerTolerance", 0);
304 ffunc->Evaluate(xc0, xc1, xc2, tolerance);
305 lowerTol = Vmath::Vsum(nq, tolerance, 1) / nq;
306 }
307
308 if (m_session->DefinesFunction("AdaptiveUpperTolerance"))
309 {
310 int nq = Exp->GetTotPoints();
311
312 // Obtain points from the the element
313 Array<OneD, NekDouble> xc0, xc1, xc2;
314 xc0 = Array<OneD, NekDouble>(nq, 0.0);
315 xc1 = Array<OneD, NekDouble>(nq, 0.0);
316 xc2 = Array<OneD, NekDouble>(nq, 0.0);
317 Exp->GetCoords(xc0, xc1, xc2);
318
319 // Evaluate function from session file
320 Array<OneD, NekDouble> tolerance(nq, 0.0);
322 m_session->GetFunction("AdaptiveUpperTolerance", 0);
323 ffunc->Evaluate(xc0, xc1, xc2, tolerance);
324 upperTol = Vmath::Vsum(nq, tolerance, 1) / nq;
325 }
326
327 // Determine what to do with the polynomial order
328 if ((error > upperTol) && (P[0] < maxP))
329 {
330 deltaP[Exp->GetGeom()->GetGlobalID()] = 1;
331 }
332 else if ((error < lowerTol) && P[0] > minP)
333 {
334 deltaP[Exp->GetGeom()->GetGlobalID()] = -1;
335 }
336 else
337 {
338 deltaP[Exp->GetGeom()->GetGlobalID()] = 0;
339 }
340 }
341
342 // Write new expansion section to the session reader and re-read graph.
343 ReplaceExpansion(fields, deltaP);
344 m_graph->ReadExpansionInfo();
345
346 // Reset GlobalLinSys Manager to avoid using too much memory
347 //
348 // @todo This could be made better by replacing individual matrices
349 // within the linear system.
352 PoolCreated(std::string("GlobalLinSys")))
353 {
356 ClearManager(std::string("GlobalLinSys"));
357 }
358
359 int chkNumber = m_equ[0]->GetCheckpointNumber();
360 int chkSteps = m_equ[0]->GetCheckpointSteps();
361
362 // Initialise driver again
364
365 // Update mapping (must be before m_equ[0]->DoInitialise();)
366 mapping->ReplaceField(m_equ[0]->UpdateFields());
367
368 // Set chkSteps to zero to avoid writing initial condition
369 m_equ[0]->SetCheckpointSteps(0);
370
371 // Initialise equation
372 m_equ[0]->DoInitialise();
373 m_equ[0]->SetInitialStep(i * numSteps);
374 m_equ[0]->SetSteps(i * numSteps + numSteps);
375 m_equ[0]->SetTime(startTime + i * period);
376 m_equ[0]->SetBoundaryConditions(startTime + i * period);
377 m_equ[0]->SetCheckpointNumber(chkNumber);
378 m_equ[0]->SetCheckpointSteps(chkSteps);
379
380 // Project solution to new expansion
381 for (int n = 0; n < nFields; n++)
382 {
383 m_equ[0]->UpdateFields()[n]->ExtractCoeffsToCoeffs(
384 fields[n], fields[n]->GetCoeffs(),
385 m_equ[0]->UpdateFields()[n]->UpdateCoeffs());
386 m_equ[0]->UpdateFields()[n]->BwdTrans(
387 m_equ[0]->UpdateFields()[n]->GetCoeffs(),
388 m_equ[0]->UpdateFields()[n]->UpdatePhys());
389 }
390
391 // Solve equation
392 m_equ[0]->DoSolve();
393 }
394
395 timer.Stop();
396
397 m_equ[0]->Output();
398
399 if (m_comm->GetRank() == 0)
400 {
401 CPUtime = timer.Elapsed().count();
402 std::cout << "-------------------------------------------" << std::endl;
403 std::cout << "Total Computation Time = " << CPUtime << "s" << std::endl;
404 std::cout << "-------------------------------------------" << std::endl;
405 }
406
407 // Evaluate and output computation time and solution accuracy.
408 // The specific format of the error output is essential for the
409 // regression tests to work.
410
411 // Evaluate L2 Error
412 for (int i = 0; i < m_equ[0]->GetNvariables(); ++i)
413 {
414 Array<OneD, NekDouble> exactsoln(m_equ[0]->GetTotPoints(), 0.0);
415
416 // Evaluate "ExactSolution" function, or zero array
417 m_equ[0]->EvaluateExactSolution(i, exactsoln, m_equ[0]->GetTime());
418
419 NekDouble vL2Error = m_equ[0]->L2Error(i, exactsoln);
420 NekDouble vLinfError = m_equ[0]->LinfError(i, exactsoln);
421
422 if (m_comm->GetRank() == 0)
423 {
424 out << "L 2 error (variable " << m_equ[0]->GetVariable(i)
425 << ") : " << vL2Error << std::endl;
426 out << "L inf error (variable " << m_equ[0]->GetVariable(i)
427 << ") : " << vLinfError << std::endl;
428 }
429 }
430}
431
432/**
433 * @brief Update EXPANSIONS tag inside XML schema to reflect new polynomial
434 * order distribution.
435 *
436 * @param fields Input fields.
437 * @param deltaP Map of polynomial order expansions
438 */
441 std::map<int, int> deltaP)
442{
443 int nExp, nDim;
444 int expdim = m_equ[0]->UpdateFields()[0]->GetGraph()->GetMeshDimension();
445
446 // Get field definitions
447 std::vector<LibUtilities::FieldDefinitionsSharedPtr> fielddefs =
448 fields[0]->GetFieldDefinitions();
449
450 if (fielddefs[0]->m_numHomogeneousDir == 1)
451 {
452 nDim = expdim + 1;
453 }
454 else
455 {
456 nDim = expdim;
457 }
458
459 // Add variables to field definition
460 for (int i = 0; i < fielddefs.size(); ++i)
461 {
462 for (int j = 0; j < fields.size(); ++j)
463 {
464 fielddefs[i]->m_fields.push_back(m_session->GetVariable(j));
465 }
466 }
467
468 // Get tinyxml objects
469 TiXmlElement *exp_tag = m_session->GetElement("NEKTAR/EXPANSIONS");
470
471 // Clear current expansions
472 exp_tag->Clear();
473
474 // Write new expansion information
475 for (int f = 0; f < fielddefs.size(); ++f)
476 {
477 nExp = fielddefs[f]->m_elementIDs.size();
478
479 //---------------------------------------------
480 // Write ELEMENTS
481 TiXmlElement *elemTag = new TiXmlElement("ELEMENTS");
482 exp_tag->LinkEndChild(elemTag);
483
484 // Write FIELDS
485 std::string fieldsString;
486 {
487 std::stringstream fieldsStringStream;
488 bool first = true;
489 for (std::vector<int>::size_type i = 0;
490 i < fielddefs[f]->m_fields.size(); i++)
491 {
492 if (!first)
493 {
494 fieldsStringStream << ",";
495 }
496 fieldsStringStream << fielddefs[f]->m_fields[i];
497 first = false;
498 }
499 fieldsString = fieldsStringStream.str();
500 }
501 elemTag->SetAttribute("FIELDS", fieldsString);
502
503 // Write SHAPE
504 std::string shapeString;
505 {
506 std::stringstream shapeStringStream;
507 shapeStringStream
508 << LibUtilities::ShapeTypeMap[fielddefs[f]->m_shapeType];
509 if (fielddefs[f]->m_numHomogeneousDir == 1)
510 {
511 shapeStringStream << "-HomogenousExp1D";
512 }
513 else if (fielddefs[f]->m_numHomogeneousDir == 2)
514 {
515 shapeStringStream << "-HomogenousExp2D";
516 }
517
518 shapeString = shapeStringStream.str();
519 }
520 elemTag->SetAttribute("SHAPE", shapeString);
521
522 // Write BASIS
523 std::string basisString;
524 {
525 std::stringstream basisStringStream;
526 bool first = true;
527 for (std::vector<LibUtilities::BasisType>::size_type i = 0;
528 i < fielddefs[f]->m_basis.size(); i++)
529 {
530 if (!first)
531 {
532 basisStringStream << ",";
533 }
534 basisStringStream
535 << LibUtilities::BasisTypeMap[fielddefs[f]->m_basis[i]];
536 first = false;
537 }
538 basisString = basisStringStream.str();
539 }
540 elemTag->SetAttribute("BASIS", basisString);
541
542 // Write homogeneuous length details
543 if (fielddefs[f]->m_numHomogeneousDir)
544 {
545 std::string homoLenString;
546 {
547 std::stringstream homoLenStringStream;
548 bool first = true;
549 for (int i = 0; i < fielddefs[f]->m_numHomogeneousDir; ++i)
550 {
551 if (!first)
552 {
553 homoLenStringStream << ",";
554 }
555 homoLenStringStream
556 << fielddefs[f]->m_homogeneousLengths[i];
557 first = false;
558 }
559 homoLenString = homoLenStringStream.str();
560 }
561 elemTag->SetAttribute("HOMOGENEOUSLENGTHS", homoLenString);
562 }
563
564 // Write homogeneuous planes/lines details
565 if (fielddefs[f]->m_numHomogeneousDir)
566 {
567 if (fielddefs[f]->m_homogeneousYIDs.size() > 0)
568 {
569 std::string homoYIDsString;
570 {
571 std::stringstream homoYIDsStringStream;
572 bool first = true;
573 for (int i = 0; i < fielddefs[f]->m_homogeneousYIDs.size();
574 i++)
575 {
576 if (!first)
577 {
578 homoYIDsStringStream << ",";
579 }
580 homoYIDsStringStream
581 << fielddefs[f]->m_homogeneousYIDs[i];
582 first = false;
583 }
584 homoYIDsString = homoYIDsStringStream.str();
585 }
586 elemTag->SetAttribute("HOMOGENEOUSYIDS", homoYIDsString);
587 }
588
589 if (fielddefs[f]->m_homogeneousZIDs.size() > 0)
590 {
591 std::string homoZIDsString;
592 {
593 std::stringstream homoZIDsStringStream;
594 bool first = true;
595 for (int i = 0; i < fielddefs[f]->m_homogeneousZIDs.size();
596 i++)
597 {
598 if (!first)
599 {
600 homoZIDsStringStream << ",";
601 }
602 homoZIDsStringStream
603 << fielddefs[f]->m_homogeneousZIDs[i];
604 first = false;
605 }
606 homoZIDsString = homoZIDsStringStream.str();
607 }
608 elemTag->SetAttribute("HOMOGENEOUSZIDS", homoZIDsString);
609 }
610 }
611
612 // Write NUMMODESPERDIR
613 std::string numModesString;
614 {
615 std::stringstream numModesStringStream;
616
617 numModesStringStream << "MIXORDER:";
618 bool first = true;
619 int eId;
620 Array<OneD, int> order(nDim, 0);
621 for (int n = 0; n < nExp; n++)
622 {
623 eId = fielddefs[f]->m_elementIDs[n];
624
625 for (int i = 0; i < expdim; i++)
626 {
627 order[i] = deltaP[eId];
628 }
629
630 for (int i = 0; i < nDim; i++)
631 {
632 if (fielddefs[f]->m_uniOrder)
633 {
634 order[i] += fielddefs[f]->m_numModes[i];
635 }
636 else
637 {
638 order[i] += fielddefs[f]->m_numModes[n * nDim + i];
639 }
640
641 if (!first)
642 {
643 numModesStringStream << ",";
644 }
645
646 numModesStringStream << order[i];
647 first = false;
648 }
649 }
650 numModesString = numModesStringStream.str();
651 }
652 elemTag->SetAttribute("NUMMODESPERDIR", numModesString);
653
654 // Write IDs. Should ideally look at ways of compressing this stream if
655 // just sequential.
656 elemTag->SetAttribute(
657 "ID", ParseUtils::GenerateSeqString(fielddefs[f]->m_elementIDs));
658 }
659}
660
661} // namespace Nektar::SolverUtils
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
static GLOBAL_MAPPING_EXPORT MappingSharedPtr Load(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Return a pointer to the mapping, creating it on first call.
Definition: Mapping.cpp:264
Describes the specification for a Basis.
Definition: Basis.h:45
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Defines a specification for a set of points.
Definition: Points.h:50
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
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.
A global linear system.
Definition: GlobalLinSys.h:70
static std::string GenerateSeqString(const std::vector< T > &v)
Generate a compressed comma-separated string representation of a vector of unsigned integers.
Definition: ParseUtils.h:72
static DriverSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
SOLVER_UTILS_EXPORT void ReplaceExpansion(Array< OneD, MultiRegions::ExpListSharedPtr > &fields, std::map< int, int > deltaP)
Update EXPANSIONS tag inside XML schema to reflect new polynomial order distribution.
SOLVER_UTILS_EXPORT DriverAdaptive(const LibUtilities::SessionReaderSharedPtr pSession, const SpatialDomains::MeshGraphSharedPtr pGraph)
Constructor.
SOLVER_UTILS_EXPORT void v_InitObject(std::ostream &out=std::cout) override
Virtual function for initialisation implementation.
SOLVER_UTILS_EXPORT void v_Execute(std::ostream &out=std::cout) override
Virtual function for solve implementation.
static std::string className
Name of the class.
Base class for the development of solvers.
Definition: Driver.h:65
LibUtilities::SessionReaderSharedPtr m_session
Session reader object.
Definition: Driver.h:83
virtual SOLVER_UTILS_EXPORT void v_InitObject(std::ostream &out=std::cout)
Virtual function for initialisation implementation.
Definition: Driver.cpp:82
LibUtilities::CommSharedPtr m_comm
Communication object.
Definition: Driver.h:80
SpatialDomains::MeshGraphSharedPtr m_graph
MeshGraph object.
Definition: Driver.h:89
Array< OneD, EquationSystemSharedPtr > m_equ
Equation system to solve.
Definition: Driver.h:92
Class representing a prismatic element in reference space.
Definition: StdPrismExp.h:45
GLOBAL_MAPPING_EXPORT typedef std::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
Definition: Mapping.h:57
const char *const ShapeTypeMap[SIZE_ShapeType]
Definition: ShapeType.hpp:75
const char *const BasisTypeMap[]
Definition: Foundations.hpp:44
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:125
@ P
Monomial polynomials .
Definition: BasisType.h:62
@ eOrtho_A
Principle Orthogonal Functions .
Definition: BasisType.h:42
@ eOrtho_C
Principle Orthogonal Functions .
Definition: BasisType.h:46
@ eOrtho_B
Principle Orthogonal Functions .
Definition: BasisType.h:44
std::shared_ptr< Expansion > ExpansionSharedPtr
Definition: Expansion.h:66
DriverFactory & GetDriverFactory()
Definition: Driver.cpp:64
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::shared_ptr< StdExpansion > StdExpansionSharedPtr
double NekDouble
void Vmul(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Multiply vector z = x*y.
Definition: Vmath.hpp:72
T Vsum(int n, const T *x, const int incx)
Subtract return sum(x)
Definition: Vmath.hpp:608
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.
Definition: Vmath.hpp:220
scalarT< T > abs(scalarT< T > in)
Definition: scalar.hpp:289