Nektar++
Mapping.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: Mapping.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: Abstract base class for mappings.
32//
33///////////////////////////////////////////////////////////////////////////////
34
37
38#include <boost/algorithm/string/predicate.hpp>
39
41{
42
44bool Mapping::m_init = false;
45bool Mapping::m_isDefined = false;
46
48{
49 static MappingFactory instance;
50 return instance;
51}
52
55 : m_session(pSession), m_fields(pFields)
56{
57 switch (m_fields[0]->GetExpType())
58 {
60 {
62 }
63 break;
64
66 {
68 }
69 break;
70
74 {
76 }
77 break;
78
79 default:
80 ASSERTL0(0, "Dimension not supported");
81 break;
82 }
83
85}
86
87/**
88 * This function initialises the Mapping object. It computes the coordinates
89 * and velocity coordinates, initialises the workspace variables, and calls
90 * UpdateGeomInfo, which will perform the calculations specific for each type
91 * of Mapping.
92 * @param pFields ExpList array used in the mapping
93 * @param pMapping xml element describing the mapping
94 */
96 [[maybe_unused]] const Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
97 const TiXmlElement *pMapping)
98{
99 int phystot = m_fields[0]->GetTotPoints();
100 m_fromFunction = true;
101 // Initialise variables
105 for (int i = 0; i < 3; i++)
106 {
107 m_coords[i] = Array<OneD, NekDouble>(phystot);
109 coords[i] = Array<OneD, NekDouble>(phystot);
110 }
111
112 // Check if mapping is defined as time-dependent
113 const TiXmlElement *timeDep = pMapping->FirstChildElement("TIMEDEPENDENT");
114 if (timeDep)
115 {
116 std::string sTimeDep = timeDep->GetText();
117 m_timeDependent = (boost::iequals(sTimeDep, "true")) ||
118 (boost::iequals(sTimeDep, "yes"));
119 }
120 else
121 {
122 m_timeDependent = false;
123 }
124
125 // Load coordinates
126 std::string fieldNames[3] = {"x", "y", "z"};
127 const TiXmlElement *funcNameElmt = pMapping->FirstChildElement("COORDS");
128 if (funcNameElmt)
129 {
130 m_funcName = funcNameElmt->GetText();
131 ASSERTL0(m_session->DefinesFunction(m_funcName),
132 "Function '" + m_funcName + "' not defined.");
133
134 // Get coordinates in the domain
135 m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
136
137 std::string s_FieldStr;
138 // Check if function from session file defines each component
139 // and evaluate them, otherwise use trivial transformation
140 for (int i = 0; i < 3; i++)
141 {
142 s_FieldStr = fieldNames[i];
143 if (m_session->DefinesFunction(m_funcName, s_FieldStr))
144 {
146 m_funcName);
147 if (i == 2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
148 {
149 ASSERTL0(
150 false,
151 "3DH1D does not support mapping in the z-direction.");
152 }
153 }
154 else
155 {
156 // This coordinate is not defined, so use (x^i)' = x^i
157 Vmath::Vcopy(phystot, coords[i], 1, m_coords[i], 1);
158 }
159 }
160 }
161 else
162 {
163 m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
164 for (int i = 0; i < 3; i++)
165 {
166 // Use (x^i)' = x^i as default. This can be useful if we
167 // have a time-dependent mapping, and then only the
168 // initial mapping will be trivial
169 Vmath::Vcopy(phystot, coords[i], 1, m_coords[i], 1);
170 }
171 }
172
173 // Load coordinate velocity if they are defined,
174 // otherwise use zero to make it general
175 std::string velFieldNames[3] = {"vx", "vy", "vz"};
176 const TiXmlElement *velFuncNameElmt = pMapping->FirstChildElement("VEL");
177 if (velFuncNameElmt)
178 {
179 m_velFuncName = velFuncNameElmt->GetText();
180 ASSERTL0(m_session->DefinesFunction(m_velFuncName),
181 "Function '" + m_velFuncName + "' not defined.");
182
183 std::string s_FieldStr;
184 // Check if function from session file defines each component
185 // and evaluate them, otherwise use 0
186 for (int i = 0; i < 3; i++)
187 {
188 s_FieldStr = velFieldNames[i];
189 if (m_session->DefinesFunction(m_velFuncName, s_FieldStr))
190 {
193 if (i == 2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
194 {
195 ASSERTL0(
196 false,
197 "3DH1D does not support mapping in the z-direction.");
198 }
199 }
200 else
201 {
202 // This coordinate velocity is not defined, so use 0
203 Vmath::Zero(phystot, m_coordsVel[i], 1);
204 }
205 }
206 }
207 else
208 {
209 for (int i = 0; i < 3; i++)
210 {
211 Vmath::Zero(phystot, m_coordsVel[i], 1);
212 }
213 }
214
215 // Initialise workspace variables
216 int nvel = m_nConvectiveFields;
220 for (int i = 0; i < nvel; i++)
221 {
222 m_tmp[i] = Array<OneD, NekDouble>(phystot, 0.0);
223 for (int j = 0; j < nvel; j++)
224 {
225 m_wk1[i * nvel + j] = Array<OneD, NekDouble>(phystot, 0.0);
226 m_wk2[i * nvel + j] = Array<OneD, NekDouble>(phystot, 0.0);
227 }
228 }
229
230 // Calculate information required by the particular mapping
232}
233
234/**
235 * This function replaces the expansion list contained in m_fields, and then
236 * proceeds reinitialising the Mapping with this new field.
237 * @param pFields New field to be used by the Mapping
238 */
241{
242 m_fields = pFields;
243
244 TiXmlElement *vMapping = nullptr;
245
246 if (m_session->DefinesElement("Nektar/Mapping"))
247 {
248 vMapping = m_session->GetElement("Nektar/Mapping");
250 vMapping, m_session->GetTimeLevel());
251 }
252 InitObject(pFields, vMapping);
253}
254
255/**
256 * This function is responsible for loading the Mapping, guaranteeing that a
257 * single instance of this class exists. When it is first called, it creates a
258 * Mapping and returns a pointer to it. On subsequent calls, it just returns
259 * the pointer.
260 * @param pSession Session reader object
261 * @param pFields Fields which will be used by the Mapping
262 * @return Pointer to the Mapping
263 */
267{
268 if (!m_init)
269 {
270 TiXmlElement *vMapping = nullptr;
271 std::string vType;
272 if (pSession->DefinesElement("Nektar/Mapping"))
273 {
274 vMapping = pSession->GetElement("Nektar/Mapping");
275 vType = vMapping->Attribute("TYPE");
276 m_isDefined = true;
278 vMapping, pSession->GetTimeLevel());
279 }
280 else
281 {
282 vType = "Translation";
283 }
284
286 pFields, vMapping);
287
288 m_init = true;
289 }
290
291 return m_mappingPtr;
292}
293
294/**
295 * This function should be called before writing a fld or chk file. It updates
296 * the metadata with information from the Mapping. Also, if the mapping is not
297 * defined by a function, it writes a .map file containing the coordinates and
298 * velocity of the coordinates.
299 * @param fieldMetaDataMap Metadata of the output file
300 * @param outname Name of the output file
301 */
303 const std::string &outname)
304{
305 // Only do anything if mapping exists
306 if (m_isDefined)
307 {
308 fieldMetaDataMap["MappingCartesianVel"] = std::string("False");
309 if (m_fromFunction)
310 {
311 // Add metadata
312 fieldMetaDataMap["MappingType"] = std::string("Expression");
313 fieldMetaDataMap["MappingExpression"] = m_funcName;
314 if (m_timeDependent)
315 {
316 fieldMetaDataMap["MappingVelExpression"] = m_velFuncName;
317 }
318 }
319 else
320 {
321 int expdim = m_fields[0]->GetGraph()->GetMeshDimension();
322 std::string fieldNames[3] = {"x", "y", "z"};
323 std::string velFieldNames[3] = {"vx", "vy", "vz"};
324
325 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef =
326 m_fields[0]->GetFieldDefinitions();
327 std::vector<std::vector<NekDouble>> FieldData(FieldDef.size());
328
329 int ncoeffs = m_fields[0]->GetNcoeffs();
330 Array<OneD, NekDouble> fieldcoeffs(ncoeffs);
331
332 bool wavespace = m_fields[0]->GetWaveSpace();
333 m_fields[0]->SetWaveSpace(false);
334 // copy coordinates Data into FieldData and set variable
335 for (int j = 0; j < expdim; ++j)
336 {
337 m_fields[0]->FwdTransLocalElmt(m_coords[j], fieldcoeffs);
338
339 for (int i = 0; i < FieldDef.size(); ++i)
340 {
341 // Could do a search here to find correct variable
342 FieldDef[i]->m_fields.push_back(fieldNames[j]);
343 m_fields[0]->AppendFieldData(FieldDef[i], FieldData[i],
344 fieldcoeffs);
345 }
346 }
347 if (m_timeDependent)
348 {
349 // copy coordinates velocity Data into FieldData and set
350 // variable
351 for (int j = 0; j < expdim; ++j)
352 {
353 m_fields[0]->FwdTransLocalElmt(m_coordsVel[j], fieldcoeffs);
354
355 for (int i = 0; i < FieldDef.size(); ++i)
356 {
357 // Could do a search here to find correct variable
358 FieldDef[i]->m_fields.push_back(velFieldNames[j]);
359 m_fields[0]->AppendFieldData(FieldDef[i], FieldData[i],
360 fieldcoeffs);
361 }
362 }
363 }
364
365 std::string outfile = outname;
366 outfile.erase(outfile.end() - 4, outfile.end());
367 outfile += ".map";
368
369 m_fld->Write(outfile, FieldDef, FieldData, fieldMetaDataMap);
370
371 // Write metadata to orginal output
372 fieldMetaDataMap["MappingType"] = std::string("File");
373 fieldMetaDataMap["FileName"] = outfile;
374
375 m_fields[0]->SetWaveSpace(wavespace);
376 }
377 }
378}
379
381 LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName,
382 Array<OneD, NekDouble> &pArray, const std::string &pFunctionName,
383 NekDouble pTime)
384{
385 ASSERTL0(pSession->DefinesFunction(pFunctionName),
386 "Function '" + pFunctionName + "' does not exist.");
387
389 pSession->GetFunction(pFunctionName, pFieldName);
390
391 Array<OneD, NekDouble> x0(1, 0.0);
392 Array<OneD, NekDouble> x1(1, 0.0);
393 Array<OneD, NekDouble> x2(1, 0.0);
394
395 ffunc->Evaluate(x0, x1, x2, pTime, pArray);
396}
397
400 LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName,
401 Array<OneD, NekDouble> &pArray, const std::string &pFunctionName,
402 NekDouble pTime)
403{
404 ASSERTL0(pSession->DefinesFunction(pFunctionName),
405 "Function '" + pFunctionName + "' does not exist.");
406
407 unsigned int nq = pFields[0]->GetNpoints();
408 if (pArray.size() != nq)
409 {
410 pArray = Array<OneD, NekDouble>(nq);
411 }
412
414 vType = pSession->GetFunctionType(pFunctionName, pFieldName);
416 {
420
421 pFields[0]->GetCoords(x0, x1, x2);
423 pSession->GetFunction(pFunctionName, pFieldName);
424
425 ffunc->Evaluate(x0, x1, x2, pTime, pArray);
426 }
427 else if (vType == LibUtilities::eFunctionTypeFile)
428 {
429 std::string filename =
430 pSession->GetFunctionFilename(pFunctionName, pFieldName);
431
432 std::vector<LibUtilities::FieldDefinitionsSharedPtr> FieldDef;
433 std::vector<std::vector<NekDouble>> FieldData;
434 Array<OneD, NekDouble> vCoeffs(pFields[0]->GetNcoeffs());
435 Vmath::Zero(vCoeffs.size(), vCoeffs, 1);
436
438 LibUtilities::FieldIO::CreateForFile(pSession, filename);
439 fld->Import(filename, FieldDef, FieldData);
440
441 int idx = -1;
442 for (int i = 0; i < FieldDef.size(); ++i)
443 {
444 for (int j = 0; j < FieldDef[i]->m_fields.size(); ++j)
445 {
446 if (FieldDef[i]->m_fields[j] == pFieldName)
447 {
448 idx = j;
449 }
450 }
451
452 if (idx >= 0)
453 {
454 pFields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],
455 FieldDef[i]->m_fields[idx],
456 vCoeffs);
457 }
458 else
459 {
460 std::cout << "Field " + pFieldName + " not found." << std::endl;
461 }
462 }
463 pFields[0]->BwdTrans(vCoeffs, pArray);
464 }
465}
466
467/**
468 * This function converts a contravariant vector in transformed space
469 * \f$v^{j}\f$ to the corresponding \f$\bar{v}^{i}\f$ in Cartesian (physical)
470 * space using the relation
471 * \f[\bar{v}^{i} = \frac{\partial \bar{x}^i}{\partial x^j}v^{j}\f]
472 *
473 * @param inarray Components of the vector in transformed space (\f$v^{j}\f$)
474 * @param outarray Components of the vector in Cartesian space
475 * (\f$\bar{v}^{i}\f$)
476 */
478 const Array<OneD, Array<OneD, NekDouble>> &inarray,
480{
481 if (inarray == outarray)
482 {
483 int physTot = m_fields[0]->GetTotPoints();
484 int nvel = m_nConvectiveFields;
485
486 for (int i = 0; i < nvel; i++)
487 {
488 Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
489 }
490 v_ContravarToCartesian(m_tmp, outarray);
491 }
492 else
493 {
494 v_ContravarToCartesian(inarray, outarray);
495 }
496}
497
498/**
499 * This function converts a covariant vector in transformed space
500 * \f$v_{j}\f$ to the corresponding \f$\bar{v}_{i}\f$ in Cartesian (physical)
501 * space using the relation
502 * \f[\bar{v}_{i} = \frac{\partial x^j}{\partial \bar{x}^i}v_{j}\f]
503 *
504 * @param inarray Components of the vector in transformed space (\f$v_{j}\f$)
505 * @param outarray Components of the vector in Cartesian space
506 * (\f$\bar{v}_{i}\f$)
507 */
509 const Array<OneD, Array<OneD, NekDouble>> &inarray,
511{
512 if (inarray == outarray)
513 {
514 int physTot = m_fields[0]->GetTotPoints();
515 int nvel = m_nConvectiveFields;
516
517 for (int i = 0; i < nvel; i++)
518 {
519 Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
520 }
521 v_CovarToCartesian(m_tmp, outarray);
522 }
523 else
524 {
525 v_CovarToCartesian(inarray, outarray);
526 }
527}
528
529/**
530 * This function converts a contravariant vector in Cartesian space
531 * \f$\bar{v}^{i}\f$ to the corresponding \f$v^{j}\f$ in transformed
532 * space using the relation
533 * \f[v^{j} = \frac{\partial x^j}{\partial \bar{x}^i}\bar{v}^{i}\f]
534 *
535 * @param inarray Components of the vector in Cartesian space
536 * (\f$\bar{v}^{i}\f$)
537 * @param outarray Components of the vector in transformed space (\f$v^{j}\f$)
538 */
540 const Array<OneD, Array<OneD, NekDouble>> &inarray,
542{
543 if (inarray == outarray)
544 {
545 int physTot = m_fields[0]->GetTotPoints();
546 int nvel = m_nConvectiveFields;
547
548 for (int i = 0; i < nvel; i++)
549 {
550 Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
551 }
553 }
554 else
555 {
556 v_ContravarFromCartesian(inarray, outarray);
557 }
558}
559
560/**
561 * This function converts a covariant vector in Cartesian space
562 * \f$\bar{v}_{i}\f$ to the corresponding \f$v_{j}\f$ in transformed
563 * space using the relation
564 * \f[\bar{v}_{j} = \frac{\partial \bar{x}^i}{\partial x^j}v_{i}\f]
565 *
566 * @param inarray Components of the vector in Cartesian space
567 * (\f$\bar{v}_{i}\f$)
568 * @param outarray Components of the vector in transformed space (\f$v_{j}\f$)
569 */
571 const Array<OneD, Array<OneD, NekDouble>> &inarray,
573{
574 if (inarray == outarray)
575 {
576 int physTot = m_fields[0]->GetTotPoints();
577 int nvel = m_nConvectiveFields;
578
579 for (int i = 0; i < nvel; i++)
580 {
581 Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
582 }
583 v_CovarFromCartesian(m_tmp, outarray);
584 }
585 else
586 {
587 v_CovarFromCartesian(inarray, outarray);
588 }
589}
590
591/**
592 * This function lowers the index of the contravariant vector \f$v^{i}\f$,
593 * transforming it into its associated covariant vector \f$v_{j}\f$
594 * according to the relation
595 * \f[v_{j} = g_{ij}v^{i}\f]
596 * where \f$g_{ij}\f$ is the metric tensor.
597 *
598 * @param inarray Components of the contravariant vector \f$v^{i}\f$
599 * @param outarray Components of the covariant vector \f$v_{j}\f$
600 */
603{
604 if (inarray == outarray)
605 {
606 int physTot = m_fields[0]->GetTotPoints();
607 int nvel = m_nConvectiveFields;
608
609 for (int i = 0; i < nvel; i++)
610 {
611 Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
612 }
613 v_LowerIndex(m_tmp, outarray);
614 }
615 else
616 {
617 v_LowerIndex(inarray, outarray);
618 }
619}
620
621/**
622 * This function raises the index of the covariant vector \f$v_{j}\f$,
623 * transforming it into its associated contravariant vector \f$v^{i}\f$
624 * according to the relation
625 * \f[v^{i} = g^{ij}v_{j}\f]
626 * where \f$g^{ij}\f$ is the inverse of the metric tensor.
627 *
628 * @param inarray Components of the contravariant vector \f$v^{i}\f$
629 * @param outarray Components of the covariant vector \f$v_{j}\f$
630 */
633{
634 if (inarray == outarray)
635 {
636 int physTot = m_fields[0]->GetTotPoints();
637 int nvel = m_nConvectiveFields;
638
639 for (int i = 0; i < nvel; i++)
640 {
641 Vmath::Vcopy(physTot, inarray[i], 1, m_tmp[i], 1);
642 }
643 v_RaiseIndex(m_tmp, outarray);
644 }
645 else
646 {
647 v_RaiseIndex(inarray, outarray);
648 }
649}
650
654{
655 int physTot = m_fields[0]->GetTotPoints();
656
657 out0 = Array<OneD, NekDouble>(physTot, 0.0);
658 out1 = Array<OneD, NekDouble>(physTot, 0.0);
659 out2 = Array<OneD, NekDouble>(physTot, 0.0);
660
661 Vmath::Vcopy(physTot, m_coords[0], 1, out0, 1);
662 Vmath::Vcopy(physTot, m_coords[1], 1, out1, 1);
663 Vmath::Vcopy(physTot, m_coords[2], 1, out2, 1);
664}
665
667{
668 int physTot = m_fields[0]->GetTotPoints();
669
670 for (int i = 0; i < m_nConvectiveFields; ++i)
671 {
672 outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
673 Vmath::Vcopy(physTot, m_coordsVel[i], 1, outarray[i], 1);
674 }
675}
676
678 const Array<OneD, Array<OneD, NekDouble>> &inarray,
679 Array<OneD, NekDouble> &outarray)
680{
681 int physTot = m_fields[0]->GetTotPoints();
682
683 outarray = Array<OneD, NekDouble>(physTot, 0.0);
684 if (!HasConstantJacobian())
685 {
686 // Set wavespace to false and store current value
687 bool wavespace = m_fields[0]->GetWaveSpace();
688 m_fields[0]->SetWaveSpace(false);
689
690 // Get Mapping Jacobian
691 Array<OneD, NekDouble> Jac(physTot, 0.0);
692 GetJacobian(Jac);
693
694 // Calculate inarray . grad(Jac)
695 Array<OneD, NekDouble> wk(physTot, 0.0);
696 for (int i = 0; i < m_nConvectiveFields; ++i)
697 {
698 m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], Jac, wk);
699 Vmath::Vvtvp(physTot, inarray[i], 1, wk, 1, outarray, 1, outarray,
700 1);
701 }
702 m_fields[0]->SetWaveSpace(wavespace);
703 }
704}
705
708{
709 int physTot = m_fields[0]->GetTotPoints();
710 int nvel = m_nConvectiveFields;
711
713
715
716 for (int i = 0; i < nvel; i++)
717 {
718 outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
719 for (int j = 0; j < nvel; j++)
720 {
721 Vmath::Vvtvp(physTot, g[i * nvel + j], 1, inarray[j], 1,
722 outarray[i], 1, outarray[i], 1);
723 }
724 }
725}
726
729{
730 int physTot = m_fields[0]->GetTotPoints();
731 int nvel = m_nConvectiveFields;
732
734
736
737 for (int i = 0; i < nvel; i++)
738 {
739 outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
740 for (int j = 0; j < nvel; j++)
741 {
742 Vmath::Vvtvp(physTot, g[i * nvel + j], 1, inarray[j], 1,
743 outarray[i], 1, outarray[i], 1);
744 }
745 }
746}
747
749 Array<OneD, NekDouble> &outarray)
750{
751 int physTot = m_fields[0]->GetTotPoints();
752 Array<OneD, NekDouble> wk(physTot, 0.0);
753
754 Vmath::Zero(physTot, outarray, 1);
755
756 // Set wavespace to false and store current value
757 bool wavespace = m_fields[0]->GetWaveSpace();
758 m_fields[0]->SetWaveSpace(false);
759
760 // Get Mapping Jacobian
761 Array<OneD, NekDouble> Jac(physTot, 0.0);
762 GetJacobian(Jac);
763
764 for (int i = 0; i < m_nConvectiveFields; ++i)
765 {
766 Vmath::Vmul(physTot, Jac, 1, inarray[i], 1, wk, 1); // J*Ui
767 m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[i], wk,
768 wk); // (J*Ui)_i
769 Vmath::Vadd(physTot, wk, 1, outarray, 1, outarray, 1);
770 }
771 Vmath::Vdiv(physTot, outarray, 1, Jac, 1, outarray, 1); // 1/J*(J*Ui)_i
772
773 m_fields[0]->SetWaveSpace(wavespace);
774}
775
777 const Array<OneD, Array<OneD, NekDouble>> &inarray,
778 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble alpha)
779{
780 int physTot = m_fields[0]->GetTotPoints();
781 int nvel = m_nConvectiveFields;
782
784
785 // Set wavespace to false and store current value
786 bool wavespace = m_fields[0]->GetWaveSpace();
787 m_fields[0]->SetWaveSpace(false);
788
789 // Calculate vector gradient wk2 = u^i_(,k) = du^i/dx^k + {i,jk}*u^j
791 for (int i = 0; i < nvel; i++)
792 {
793 if (nvel == 2)
794 {
795 m_fields[0]->PhysDeriv(inarray[i], m_wk2[i * nvel + 0],
796 m_wk2[i * nvel + 1]);
797 }
798 else
799 {
800 m_fields[0]->PhysDeriv(inarray[i], m_wk2[i * nvel + 0],
801 m_wk2[i * nvel + 1], m_wk2[i * nvel + 2]);
802 }
803 for (int k = 0; k < nvel; k++)
804 {
805 Vmath::Vadd(physTot, m_wk1[i * nvel + k], 1, m_wk2[i * nvel + k], 1,
806 m_wk1[i * nvel + k], 1);
807 }
808 }
809 // Calculate wk1 = A^(ij) = g^(jk)*u^i_(,k)
810 for (int i = 0; i < nvel; i++)
811 {
812 for (int k = 0; k < nvel; k++)
813 {
814 tmp[k] = m_wk1[i * nvel + k];
815 }
816 RaiseIndex(tmp, m_tmp);
817 for (int j = 0; j < nvel; j++)
818 {
819 Vmath::Vcopy(physTot, m_tmp[j], 1, m_wk1[i * nvel + j], 1);
820 }
821 }
822 //
823 // Calculate L(U)^i = (A^(ij))_(,j) - alpha*d^2(u^i)/dx^jdx^j
824 //
825
826 // Step 1 :
827 // d(A^(ij) - alpha*du^i/dx^j)/d(x^j)
828 for (int i = 0; i < nvel; i++)
829 {
830 outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
831 for (int j = 0; j < nvel; j++)
832 {
833 Vmath::Smul(physTot, alpha, m_wk2[i * nvel + j], 1, m_tmp[0], 1);
834 Vmath::Vsub(physTot, m_wk1[i * nvel + j], 1, m_tmp[0], 1, m_tmp[0],
835 1);
836
837 m_fields[0]->PhysDeriv(MultiRegions::DirCartesianMap[j], m_tmp[0],
838 m_tmp[1]);
839 Vmath::Vadd(physTot, outarray[i], 1, m_tmp[1], 1, outarray[i], 1);
840 }
841 }
842
843 // Step 2: d(A^(ij))/d(x^j) + {j,pj}*A^(ip)
844 for (int i = 0; i < nvel; i++)
845 {
846 for (int p = 0; p < nvel; p++)
847 {
848 tmp[p] = m_wk1[i * nvel + p];
849 }
851 for (int j = 0; j < nvel; j++)
852 {
853 Vmath::Vadd(physTot, outarray[i], 1, m_wk2[j * nvel + j], 1,
854 outarray[i], 1);
855 }
856 }
857
858 // Step 3: d(A^(ij))/d(x^j) + {j,pj}*A^(ip) + {i,pj} A^(pj)
859 for (int j = 0; j < nvel; j++)
860 {
861 for (int p = 0; p < nvel; p++)
862 {
863 tmp[p] = m_wk1[p * nvel + j];
864 }
866 for (int i = 0; i < nvel; i++)
867 {
868 Vmath::Vadd(physTot, outarray[i], 1, m_wk2[i * nvel + j], 1,
869 outarray[i], 1);
870 }
871 }
872
873 // Restore value of wavespace
874 m_fields[0]->SetWaveSpace(wavespace);
875}
876
879{
880 int physTot = m_fields[0]->GetTotPoints();
881 int nvel = m_nConvectiveFields;
882
883 // Declare variables
884 outarray = Array<OneD, Array<OneD, NekDouble>>(nvel * nvel * nvel);
886 for (int i = 0; i < nvel * nvel * nvel; i++)
887 {
888 outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
889 }
890
891 // Set wavespace to false and store current value
892 bool wavespace = m_fields[0]->GetWaveSpace();
893 m_fields[0]->SetWaveSpace(false);
894
895 // Calculate vector gradient u^i_(,j) = du^i/dx^j + {i,pj}*u^p
897 for (int i = 0; i < nvel; i++)
898 {
899 if (nvel == 2)
900 {
901 m_fields[0]->PhysDeriv(inarray[i], m_wk2[i * nvel + 0],
902 m_wk2[i * nvel + 1]);
903 }
904 else
905 {
906 m_fields[0]->PhysDeriv(inarray[i], m_wk2[i * nvel + 0],
907 m_wk2[i * nvel + 1], m_wk2[i * nvel + 2]);
908 }
909 for (int j = 0; j < nvel; j++)
910 {
911 Vmath::Vadd(physTot, m_wk1[i * nvel + j], 1, m_wk2[i * nvel + j], 1,
912 m_wk1[i * nvel + j], 1);
913 }
914 }
915
916 //
917 // Calculate (u^i_,j),k
918 //
919
920 // Step 1 : d(u^i_,j))/d(x^k)
921 for (int i = 0; i < nvel; i++)
922 {
923 for (int j = 0; j < nvel; j++)
924 {
925 if (nvel == 2)
926 {
927 m_fields[0]->PhysDeriv(
928 m_wk1[i * nvel + j],
929 outarray[i * nvel * nvel + j * nvel + 0],
930 outarray[i * nvel * nvel + j * nvel + 1]);
931 }
932 else
933 {
934 m_fields[0]->PhysDeriv(
935 m_wk1[i * nvel + j],
936 outarray[i * nvel * nvel + j * nvel + 0],
937 outarray[i * nvel * nvel + j * nvel + 1],
938 outarray[i * nvel * nvel + j * nvel + 2]);
939 }
940 }
941 }
942
943 // Step 2: d(u^i_,j)/d(x^k) - {p,jk}*u^i_,p
944 for (int i = 0; i < nvel; i++)
945 {
946 for (int p = 0; p < nvel; p++)
947 {
948 tmp[p] = m_wk1[i * nvel + p];
949 }
951 for (int j = 0; j < nvel; j++)
952 {
953 for (int k = 0; k < nvel; k++)
954 {
955 Vmath::Vsub(physTot, outarray[i * nvel * nvel + j * nvel + k],
956 1, m_wk2[j * nvel + k], 1,
957 outarray[i * nvel * nvel + j * nvel + k], 1);
958 }
959 }
960 }
961
962 // Step 3: d(u^i_,j)/d(x^k) - {p,jk}*u^i_,p + {i,pk} u^p_,j
963 for (int j = 0; j < nvel; j++)
964 {
965 for (int p = 0; p < nvel; p++)
966 {
967 tmp[p] = m_wk1[p * nvel + j];
968 }
970 for (int i = 0; i < nvel; i++)
971 {
972 for (int k = 0; k < nvel; k++)
973 {
974 Vmath::Vadd(physTot, outarray[i * nvel * nvel + j * nvel + k],
975 1, m_wk2[i * nvel + k], 1,
976 outarray[i * nvel * nvel + j * nvel + k], 1);
977 }
978 }
979 }
980
981 // Restore value of wavespace
982 m_fields[0]->SetWaveSpace(wavespace);
983}
984
987 const bool generalized)
988{
989 int physTot = m_fields[0]->GetTotPoints();
990 int nvel = m_nConvectiveFields;
991
992 // Set wavespace to false and store current value
993 bool wavespace = m_fields[0]->GetWaveSpace();
994 m_fields[0]->SetWaveSpace(false);
995
996 // For implicit treatment of viscous terms, we want the generalized
997 // curlcurl and for explicit treatment, we want the cartesian one.
998 if (generalized)
999 {
1000 // Get the second derivatives u^i_{,jk}
1002 Array<OneD, Array<OneD, NekDouble>> ddU(nvel * nvel * nvel);
1003 gradgradU(inarray, ddU);
1004
1005 // Raise index to obtain A^{ip}_{k} = g^pj u^i_{,jk}
1006 for (int i = 0; i < nvel; ++i)
1007 {
1008 for (int k = 0; k < nvel; ++k)
1009 {
1010 // Copy to wk
1011 for (int j = 0; j < nvel; ++j)
1012 {
1013 tmp[j] = ddU[i * nvel * nvel + j * nvel + k];
1014 }
1015 RaiseIndex(tmp, m_tmp);
1016 for (int p = 0; p < nvel; ++p)
1017 {
1018 Vmath::Vcopy(physTot, m_tmp[p], 1,
1019 ddU[i * nvel * nvel + p * nvel + k], 1);
1020 }
1021 }
1022 }
1023 // The curlcurl is g^ji u^k_{kj} - g^jk u^i_kj = A^{ki}_k - A^{ik}_k
1024 for (int i = 0; i < nvel; ++i)
1025 {
1026 outarray[i] = Array<OneD, NekDouble>(physTot, 0.0);
1027 for (int k = 0; k < nvel; ++k)
1028 {
1029 Vmath::Vadd(physTot, outarray[i], 1,
1030 ddU[k * nvel * nvel + i * nvel + k], 1, outarray[i],
1031 1);
1032 Vmath::Vsub(physTot, outarray[i], 1,
1033 ddU[i * nvel * nvel + k * nvel + k], 1, outarray[i],
1034 1);
1035 }
1036 }
1037 }
1038 else
1039 {
1040 m_fields[0]->CurlCurl(inarray, outarray);
1041 }
1042
1043 // Restore value of wavespace
1044 m_fields[0]->SetWaveSpace(wavespace);
1045}
1046
1048{
1049 int physTot = m_fields[0]->GetTotPoints();
1050 int nvel = m_nConvectiveFields;
1051 int nfields = m_fields.size();
1052 int nbnds = m_fields[0]->GetBndConditions().size();
1053
1054 // Declare variables
1057
1058 Array<OneD, bool> isDirichlet(nfields);
1059
1060 Array<OneD, Array<OneD, NekDouble>> values(nfields);
1061 for (int i = 0; i < nfields; i++)
1062 {
1063 values[i] = Array<OneD, NekDouble>(physTot, 0.0);
1064 }
1065
1068 Array<OneD, Array<OneD, NekDouble>> coordVel(nvel);
1069 for (int i = 0; i < nvel; i++)
1070 {
1071 tmp[i] = Array<OneD, NekDouble>(physTot, 0.0);
1072 tmp2[i] = Array<OneD, NekDouble>(physTot, 0.0);
1073 coordVel[i] = Array<OneD, NekDouble>(physTot, 0.0);
1074 }
1075
1076 // Get coordinates velocity in transformed system (for MovingBody regions)
1077 GetCoordVelocity(tmp);
1078 ContravarFromCartesian(tmp, coordVel);
1079
1080 // Get Cartesian coordinates for evaluating boundary conditions
1082 for (int dir = 0; dir < 3; dir++)
1083 {
1084 coords[dir] = Array<OneD, NekDouble>(physTot, 0.0);
1085 }
1086 GetCartesianCoordinates(coords[0], coords[1], coords[2]);
1087
1088 // Loop boundary conditions looking for Dirichlet bc's
1089 for (int n = 0; n < nbnds; ++n)
1090 {
1091 // Evaluate original Dirichlet boundary conditions in whole domain
1092 for (int i = 0; i < nfields; ++i)
1093 {
1094 BndConds = m_fields[i]->GetBndConditions();
1095 BndExp = m_fields[i]->GetBndCondExpansions();
1096 if (BndConds[n]->GetBoundaryConditionType() ==
1098 {
1099 isDirichlet[i] = true;
1100 // If we have the a velocity component
1101 // check if all vel bc's are also Dirichlet
1102 if (i < nvel)
1103 {
1104 for (int j = 0; j < nvel; ++j)
1105 {
1107 ->GetBndConditions()[n]
1108 ->GetBoundaryConditionType() ==
1110 "Mapping only supported when all velocity "
1111 "components have the same type of boundary "
1112 "conditions");
1113 }
1114 }
1115 // Check if bc is time-dependent
1116 ASSERTL0(!BndConds[n]->IsTimeDependent(),
1117 "Time-dependent Dirichlet boundary conditions not "
1118 "supported with mapping yet.");
1119
1120 // Get boundary condition
1121 LibUtilities::Equation condition =
1122 std::static_pointer_cast<
1124 ->m_dirichletCondition;
1125 // Evaluate
1126 condition.Evaluate(coords[0], coords[1], coords[2], time,
1127 values[i]);
1128 }
1129 else
1130 {
1131 isDirichlet[i] = false;
1132 }
1133 }
1134 // Convert velocity vector to transformed system
1135 if (isDirichlet[0])
1136 {
1137 for (int i = 0; i < nvel; ++i)
1138 {
1139 Vmath::Vcopy(physTot, values[i], 1, tmp[i], 1);
1140 }
1141 ContravarFromCartesian(tmp, tmp2);
1142 for (int i = 0; i < nvel; ++i)
1143 {
1144 Vmath::Vcopy(physTot, tmp2[i], 1, values[i], 1);
1145 }
1146 }
1147
1148 // Now, project result to boundary
1149 for (int i = 0; i < nfields; ++i)
1150 {
1151 BndConds = m_fields[i]->GetBndConditions();
1152 BndExp = m_fields[i]->GetBndCondExpansions();
1153
1154 if (BndConds[n]->GetUserDefined() == "" ||
1155 BndConds[n]->GetUserDefined() == "MovingBody")
1156 {
1157 m_fields[i]->ExtractPhysToBnd(n, values[i],
1158 BndExp[n]->UpdatePhys());
1159
1160 // Apply MovingBody correction
1161 if ((i < nvel) && BndConds[n]->GetUserDefined() == "MovingBody")
1162 {
1163 // Get coordinate velocity on boundary
1164 Array<OneD, NekDouble> coordVelBnd(
1165 BndExp[n]->GetTotPoints());
1166 m_fields[i]->ExtractPhysToBnd(n, coordVel[i], coordVelBnd);
1167
1168 // Apply correction
1169 Vmath::Vadd(BndExp[n]->GetTotPoints(), coordVelBnd, 1,
1170 BndExp[n]->UpdatePhys(), 1,
1171 BndExp[n]->UpdatePhys(), 1);
1172 }
1173 }
1174 }
1175 }
1176
1177 // Finally, perform FwdTrans in all fields
1178 for (int i = 0; i < m_fields.size(); ++i)
1179 {
1180 // Get boundary condition information
1181 BndConds = m_fields[i]->GetBndConditions();
1182 BndExp = m_fields[i]->GetBndCondExpansions();
1183
1184 for (int n = 0; n < BndConds.size(); ++n)
1185 {
1186 if (BndConds[n]->GetBoundaryConditionType() ==
1188 {
1189 if (BndConds[n]->GetUserDefined() == "" ||
1190 BndConds[n]->GetUserDefined() == "MovingBody")
1191 {
1192 if (m_fields[i]->GetExpType() == MultiRegions::e3DH1D)
1193 {
1194 BndExp[n]->SetWaveSpace(false);
1195 }
1196
1197 BndExp[n]->FwdTransBndConstrained(
1198 BndExp[n]->GetPhys(), BndExp[n]->UpdateCoeffs());
1199 }
1200 }
1201 }
1202 }
1203}
1204
1206 const NekDouble time, const Array<OneD, Array<OneD, NekDouble>> &coords,
1207 const Array<OneD, Array<OneD, NekDouble>> &coordsVel)
1208{
1209 if (m_fromFunction)
1210 {
1211 std::string s_FieldStr;
1212 std::string fieldNames[3] = {"x", "y", "z"};
1213 std::string velFieldNames[3] = {"vx", "vy", "vz"};
1214 // Check if function from session file defines each component
1215 // and evaluate them, otherwise there is no need to update
1216 // coords
1217 for (int i = 0; i < 3; i++)
1218 {
1219 s_FieldStr = fieldNames[i];
1220 if (m_session->DefinesFunction(m_funcName, s_FieldStr))
1221 {
1223 m_funcName, time);
1224 if (i == 2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
1225 {
1226 ASSERTL0(
1227 false,
1228 "3DH1D does not support mapping in the z-direction.");
1229 }
1230 }
1231 s_FieldStr = velFieldNames[i];
1232 if (m_session->DefinesFunction(m_velFuncName, s_FieldStr))
1233 {
1234 EvaluateFunction(m_fields, m_session, s_FieldStr,
1235 m_coordsVel[i], m_velFuncName, time);
1236 if (i == 2 && m_fields[0]->GetExpType() == MultiRegions::e3DH1D)
1237 {
1238 ASSERTL0(
1239 false,
1240 "3DH1D does not support mapping in the z-direction.");
1241 }
1242 }
1243 }
1244 }
1245 else
1246 {
1247 int physTot = m_fields[0]->GetTotPoints();
1248 int nvel = m_nConvectiveFields;
1249 // Copy coordinates
1250 for (int i = 0; i < 3; i++)
1251 {
1252 Vmath::Vcopy(physTot, coords[i], 1, m_coords[i], 1);
1253 }
1254
1255 for (int i = 0; i < nvel; i++)
1256 {
1257 Vmath::Vcopy(physTot, coordsVel[i], 1, m_coordsVel[i], 1);
1258 }
1259 }
1260
1261 // Update the information required by the specific mapping
1263}
1264
1265} // namespace Nektar::GlobalMapping
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
bool m_fromFunction
Flag defining if the Mapping is defined by a function.
Definition: Mapping.h:431
int m_nConvectiveFields
Number of velocity components.
Definition: Mapping.h:418
static MappingSharedPtr m_mappingPtr
Definition: Mapping.h:434
virtual GLOBAL_MAPPING_EXPORT void v_GetCartesianCoordinates(Array< OneD, NekDouble > &out0, Array< OneD, NekDouble > &out1, Array< OneD, NekDouble > &out2)
Definition: Mapping.cpp:651
GLOBAL_MAPPING_EXPORT void LowerIndex(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Lower index of vector: .
Definition: Mapping.cpp:601
GLOBAL_MAPPING_EXPORT void EvaluateTimeFunction(LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
Definition: Mapping.cpp:380
LibUtilities::SessionReaderSharedPtr m_session
Session reader.
Definition: Mapping.h:406
GLOBAL_MAPPING_EXPORT void ContravarToCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Convert a contravariant vector to the Cartesian system.
Definition: Mapping.cpp:477
GLOBAL_MAPPING_EXPORT void GetInvMetricTensor(Array< OneD, Array< OneD, NekDouble > > &outarray)
Get the inverse of metric tensor .
Definition: Mapping.h:188
virtual GLOBAL_MAPPING_EXPORT void v_UpdateBCs(const NekDouble time)
Definition: Mapping.cpp:1047
GLOBAL_MAPPING_EXPORT void GetCartesianCoordinates(Array< OneD, NekDouble > &out0, Array< OneD, NekDouble > &out1, Array< OneD, NekDouble > &out2)
Get the Cartesian coordinates in the field.
Definition: Mapping.h:139
GLOBAL_MAPPING_EXPORT void gradgradU(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Second order covariant derivatives of a contravariant vector.
Definition: Mapping.h:310
GLOBAL_MAPPING_EXPORT bool IsTimeDependent()
Get flag defining if mapping is time-dependent.
Definition: Mapping.h:342
Array< OneD, Array< OneD, NekDouble > > m_wk2
Definition: Mapping.h:440
virtual GLOBAL_MAPPING_EXPORT void v_CurlCurlField(Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const bool generalized)
Definition: Mapping.cpp:985
GLOBAL_MAPPING_EXPORT void UpdateGeomInfo()
Recompute the metric terms of the Mapping.
Definition: Mapping.h:399
Array< OneD, Array< OneD, NekDouble > > m_coords
Array with the Cartesian coordinates.
Definition: Mapping.h:412
Array< OneD, Array< OneD, NekDouble > > m_coordsVel
Array with the velocity of the coordinates.
Definition: Mapping.h:414
virtual GLOBAL_MAPPING_EXPORT void v_RaiseIndex(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: Mapping.cpp:727
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Definition: Mapping.h:410
std::string m_velFuncName
Name of the function containing the velocity of the coordinates.
Definition: Mapping.h:423
GLOBAL_MAPPING_EXPORT void InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pMapping)
Initialise the mapping object.
Definition: Mapping.h:81
virtual GLOBAL_MAPPING_EXPORT void v_LowerIndex(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: Mapping.cpp:706
GLOBAL_MAPPING_EXPORT void Output(LibUtilities::FieldMetaDataMap &fieldMetaDataMap, const std::string &outname)
Output function called when a chk or fld file is written.
Definition: Mapping.cpp:302
GLOBAL_MAPPING_EXPORT void GetMetricTensor(Array< OneD, Array< OneD, NekDouble > > &outarray)
Get the metric tensor .
Definition: Mapping.h:181
Array< OneD, Array< OneD, NekDouble > > m_wk1
Definition: Mapping.h:439
GLOBAL_MAPPING_EXPORT void ContravarFromCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Convert a contravariant vector to the transformed system.
Definition: Mapping.cpp:539
virtual GLOBAL_MAPPING_EXPORT void v_GetCoordVelocity(Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: Mapping.cpp:666
GLOBAL_MAPPING_EXPORT void EvaluateFunction(Array< OneD, MultiRegions::ExpListSharedPtr > pFields, LibUtilities::SessionReaderSharedPtr pSession, std::string pFieldName, Array< OneD, NekDouble > &pArray, const std::string &pFunctionName, NekDouble pTime=NekDouble(0))
Definition: Mapping.cpp:398
virtual GLOBAL_MAPPING_EXPORT void v_gradgradU(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Definition: Mapping.cpp:877
virtual GLOBAL_MAPPING_EXPORT void v_ContravarFromCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)=0
LibUtilities::FieldIOSharedPtr m_fld
Definition: Mapping.h:408
virtual GLOBAL_MAPPING_EXPORT void v_Divergence(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition: Mapping.cpp:748
virtual GLOBAL_MAPPING_EXPORT void v_CovarFromCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)=0
GLOBAL_MAPPING_EXPORT bool HasConstantJacobian()
Get flag defining if mapping has constant Jacobian.
Definition: Mapping.h:366
GLOBAL_MAPPING_EXPORT void CovarToCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Convert a covariant vector to the Cartesian system.
Definition: Mapping.cpp:508
GLOBAL_MAPPING_EXPORT void ReplaceField(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Replace the Expansion List used by the mapping.
Definition: Mapping.cpp:239
virtual GLOBAL_MAPPING_EXPORT void v_CovarToCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)=0
GLOBAL_MAPPING_EXPORT void RaiseIndex(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Raise index of vector: .
Definition: Mapping.cpp:631
GLOBAL_MAPPING_EXPORT void ApplyChristoffelContravar(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Apply the Christoffel symbols to a contravariant vector.
Definition: Mapping.h:215
virtual GLOBAL_MAPPING_EXPORT void v_DotGradJacobian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, NekDouble > &outarray)
Definition: Mapping.cpp:677
GLOBAL_MAPPING_EXPORT void GetCoordVelocity(Array< OneD, Array< OneD, NekDouble > > &outarray)
Obtain the velocity of the coordinates.
Definition: Mapping.h:248
std::string m_funcName
Name of the function containing the coordinates.
Definition: Mapping.h:421
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
GLOBAL_MAPPING_EXPORT void CovarFromCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Convert a covariant vector to the transformed system.
Definition: Mapping.cpp:570
GLOBAL_MAPPING_EXPORT void GetJacobian(Array< OneD, NekDouble > &outarray)
Get the Jacobian of the transformation.
Definition: Mapping.h:159
GLOBAL_MAPPING_EXPORT Mapping(const LibUtilities::SessionReaderSharedPtr &pSession, const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields)
Constructor.
Definition: Mapping.cpp:53
virtual GLOBAL_MAPPING_EXPORT void v_VelocityLaplacian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble alpha)
Definition: Mapping.cpp:776
virtual GLOBAL_MAPPING_EXPORT void v_ContravarToCartesian(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)=0
virtual GLOBAL_MAPPING_EXPORT void v_InitObject(const Array< OneD, MultiRegions::ExpListSharedPtr > &pFields, const TiXmlElement *pMapping)
Definition: Mapping.cpp:95
GLOBAL_MAPPING_EXPORT void ApplyChristoffelCovar(const Array< OneD, Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray)
Apply the Christoffel symbols to a covariant vector.
Definition: Mapping.h:233
virtual GLOBAL_MAPPING_EXPORT void v_UpdateMapping(const NekDouble time, const Array< OneD, Array< OneD, NekDouble > > &coords=NullNekDoubleArrayOfArray, const Array< OneD, Array< OneD, NekDouble > > &coordsVel=NullNekDoubleArrayOfArray)
Definition: Mapping.cpp:1205
Array< OneD, Array< OneD, NekDouble > > m_tmp
Definition: Mapping.h:441
bool m_timeDependent
Flag defining if the Mapping is time-dependent.
Definition: Mapping.h:429
static std::shared_ptr< FieldIO > CreateForFile(const LibUtilities::SessionReaderSharedPtr session, const std::string &filename)
Construct a FieldIO object for a given input filename.
Definition: FieldIO.cpp:223
static std::shared_ptr< FieldIO > CreateDefault(const LibUtilities::SessionReaderSharedPtr session)
Returns an object for the default FieldIO method.
Definition: FieldIO.cpp:194
Provides a generic Factory class.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
static void GetXMLElementTimeLevel(TiXmlElement *&element, const size_t timeLevel, const bool enableCheck=true)
Get XML elment time level (Parallel-in-Time)
GLOBAL_MAPPING_EXPORT typedef std::shared_ptr< Mapping > MappingSharedPtr
A shared pointer to a Mapping object.
Definition: Mapping.h:57
MappingFactory & GetMappingFactory()
Declaration of the mapping factory singleton.
Definition: Mapping.cpp:47
std::shared_ptr< FieldIO > FieldIOSharedPtr
Definition: FieldIO.h:322
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:50
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< Equation > EquationSharedPtr
Definition: Equation.h:125
MultiRegions::Direction const DirCartesianMap[]
Definition: ExpList.h:87
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
void Vvtvp(int n, const T *w, const int incw, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
vvtvp (vector times vector plus vector): z = w*x + y
Definition: Vmath.hpp:366
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.
Definition: Vmath.hpp:180
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*x.
Definition: Vmath.hpp:100
void Vdiv(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:126
void Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825
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