Nektar++
ImageWarpingSystem.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2//
3// File: ImageWarpingSystem.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: Image warping solve routines
32//
33///////////////////////////////////////////////////////////////////////////////
34
37
38namespace Nektar
39{
40
43 "ImageWarpingSystem", ImageWarpingSystem::create,
44 "Image warping system.");
45
46/**
47 *
48 */
52 : UnsteadySystem(pSession, pGraph), AdvectionSystem(pSession, pGraph)
53{
54}
55
56/**
57 *
58 */
59void ImageWarpingSystem::v_InitObject(bool DeclareField)
60{
62
63 // Define Velocity fields
65 int nq = m_fields[0]->GetNpoints();
66
67 for (int i = 0; i < m_spacedim; ++i)
68 {
70 }
71
72 // Bit of a hack: redefine u/v fields so they are continuous for
73 // Helmholtz solve.
76 m_session, m_graph, m_session->GetVariable(2));
77 m_fields[2] = fld;
79 *fld, m_graph, m_session->GetVariable(3));
80
81 // Tell UnsteadySystem to only integrate first two fields (i.e. I and
82 // phi).
83 m_intVariables.push_back(0);
84 m_intVariables.push_back(1);
85
86 // Define the normal velocity fields
87 if (m_fields[0]->GetTrace())
88 {
90 }
91
92 std::string advName;
93 std::string riemName;
94 m_session->LoadSolverInfo("AdvectionType", advName, "WeakDG");
98 m_session->LoadSolverInfo("UpwindType", riemName, "Upwind");
100 riemName, m_session);
102 this);
103
104 m_advObject->SetRiemannSolver(m_riemannSolver);
105 m_advObject->InitObject(m_session, m_fields);
106
108 {
111 }
112 else
113 {
114 ASSERTL0(false, "Implicit unsteady Advection not set up.");
115 }
116}
117
118/**
119 *
120 */
122 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
124 [[maybe_unused]] const NekDouble time)
125{
126 int npoints = GetNpoints();
127 int ncoeffs = inarray[0].size();
129
130 // Load parameter alpha.
131 m_session->LoadParameter("Alpha", m_alpha);
132
134 "CG not implemented yet.");
135
136 // Set up storage arrays.
137 Array<OneD, NekDouble> tmp(npoints);
138 Array<OneD, NekDouble> alloc(3 * npoints);
139 Array<OneD, NekDouble> dIdx1(alloc);
140 Array<OneD, NekDouble> dIdx2(alloc + npoints);
141 Array<OneD, NekDouble> dIdx3(alloc + 2 * npoints);
142
143 // Calculate grad I.
144 m_fields[0]->PhysDeriv(inarray[0], dIdx1, dIdx2);
145
146 // Set factors.
147 // TODO: Check - should be -1?
149
150 // Multiply by phi, and perform Helmholtz solve to calculate the
151 // advection velocity field.
152 for (int i = 0; i < 2; ++i)
153 {
154 Vmath::Vmul(npoints, &alloc[i * npoints], 1, inarray[1].data(), 1,
155 m_fields[i + 2]->UpdatePhys().data(), 1);
156 Vmath::Smul(npoints, 1 / m_alpha / m_alpha,
157 m_fields[i + 2]->GetPhys().data(), 1,
158 m_fields[i + 2]->UpdatePhys().data(), 1);
159 m_fields[i + 2]->HelmSolve(m_fields[i + 2]->GetPhys(),
160 m_fields[i + 2]->UpdateCoeffs(), factors);
161 m_fields[i + 2]->BwdTrans(m_fields[i + 2]->GetCoeffs(), m_velocity[i]);
162 }
163
164 // Calculate the weak advection operator for I and phi - result is put
165 // in WeakAdv and is in physical space.
166 m_advObject->Advect(2, m_fields, m_velocity, inarray, outarray, 0.0);
167 for (int i = 0; i < 2; ++i)
168 {
169 Vmath::Neg(npoints, outarray[i], 1);
170 }
171
172 // Calculate du/dx -> dIdx1, dv/dy -> dIdx2.
173 m_fields[2]->PhysDeriv(m_velocity[0], dIdx1, dIdx3);
174 m_fields[3]->PhysDeriv(m_velocity[1], dIdx3, dIdx2);
175
176 // Calculate RHS = I*div(u) = I*du/dx + I*dv/dy -> dIdx1.
177 Vmath::Vvtvvtp(npoints, dIdx1.data(), 1, inarray[0].data(), 1, dIdx2.data(),
178 1, inarray[0].data(), 1, dIdx1.data(), 1);
179
180 // Take inner product to get to coefficient space.
181 Array<OneD, NekDouble> tmp2(ncoeffs);
182 m_fields[0]->IProductWRTBase(dIdx1, tmp2);
183
184 // Multiply by elemental inverse mass matrix, backwards transform
185 m_fields[0]->MultiplyByElmtInvMass(tmp2, tmp2);
186 m_fields[0]->BwdTrans(tmp2, tmp);
187 Vmath::Vadd(npoints, outarray[0], 1, tmp, 1, outarray[0], 1);
188}
189
190/**
191 *
192 */
194 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
195 Array<OneD, Array<OneD, NekDouble>> &outarray, const NekDouble time)
196{
197 int nvariables = inarray.size();
199
200 switch (m_projectionType)
201 {
203 {
204 // Just copy over array
205 if (inarray != outarray)
206 {
207 int npoints = GetNpoints();
208
209 for (int i = 0; i < nvariables; ++i)
210 {
211 Vmath::Vcopy(npoints, inarray[i], 1, outarray[i], 1);
212 }
213 }
214 }
215 break;
216 default:
217 ASSERTL0(false, "Unknown projection scheme");
218 break;
219 }
220}
221
222/**
223 * @brief Get the normal velocity
224 */
226{
227 // Number of trace (interface) points
228 int nTracePts = GetTraceNpoints();
229
230 // Auxiliary variable to compute the normal velocity
231 Array<OneD, NekDouble> tmp(nTracePts);
232
233 // Reset the normal velocity
234 Vmath::Zero(nTracePts, m_traceVn, 1);
235
236 for (int i = 0; i < m_velocity.size(); ++i)
237 {
238 m_fields[0]->ExtractTracePhys(m_velocity[i], tmp);
239
240 Vmath::Vvtvp(nTracePts, m_traceNormals[i], 1, tmp, 1, m_traceVn, 1,
241 m_traceVn, 1);
242 }
243
244 return m_traceVn;
245}
246
247/**
248 *
249 */
251 const Array<OneD, Array<OneD, NekDouble>> &physfield,
253{
254 ASSERTL1(flux[0].size() == m_velocity.size(),
255 "Dimension of flux array and velocity array do not match");
256
257 int nq = physfield[0].size();
258
259 for (int i = 0; i < flux.size(); ++i)
260 {
261 for (int j = 0; j < flux[0].size(); ++j)
262 {
263 Vmath::Vmul(nq, physfield[i], 1, m_velocity[j], 1, flux[i][j], 1);
264 }
265 }
266}
267
268/**
269 *
270 */
272{
274}
275} // namespace Nektar
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:208
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode....
Definition: ErrorUtil.hpp:242
ImageWarpingSystem(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
SolverUtils::RiemannSolverSharedPtr m_riemannSolver
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const SpatialDomains::MeshGraphSharedPtr &pGraph)
Creates an instance of this class.
Array< OneD, Array< OneD, NekDouble > > m_velocity
void GetFluxVector(const Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, Array< OneD, NekDouble > > > &flux)
void v_GenerateSummary(SolverUtils::SummaryList &s) override
Print a summary of time stepping parameters.
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
static std::string className
Name of class.
Array< OneD, NekDouble > & GetNormalVelocity()
Get the normal velocity.
void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
Array< OneD, NekDouble > m_traceVn
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
tBaseSharedPtr CreateInstance(tKey idKey, tParam... args)
Create an instance of the class referred to by idKey.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
static std::shared_ptr< DataType > AllocateSharedPtr(const Args &...args)
Allocate a shared pointer from the memory pool.
A base class for PDEs which include an advection component.
SolverUtils::AdvectionSharedPtr m_advObject
Advection term.
SOLVER_UTILS_EXPORT void v_InitObject(bool DeclareField=true) override
Initialisation object for EquationSystem.
int m_spacedim
Spatial dimension (>= expansion dim).
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
SOLVER_UTILS_EXPORT int GetNpoints()
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
Base class for unsteady solvers.
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s) override
Print a summary of time stepping parameters.
std::shared_ptr< SessionReader > SessionReaderSharedPtr
std::shared_ptr< ContField > ContFieldSharedPtr
Definition: ContField.h:268
AdvectionFactory & GetAdvectionFactory()
Gets the factory for initialising advection objects.
Definition: Advection.cpp:43
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
EquationSystemFactory & GetEquationSystemFactory()
RiemannSolverFactory & GetRiemannSolverFactory()
std::shared_ptr< MeshGraph > MeshGraphSharedPtr
Definition: MeshGraph.h:174
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:430
StdRegions::ConstFactorMap factors
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 Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.hpp:292
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 Zero(int n, T *x, const int incx)
Zero vector.
Definition: Vmath.hpp:273
void Vvtvvtp(int n, const T *v, int incv, const T *w, int incw, const T *x, int incx, const T *y, int incy, T *z, int incz)
vvtvvtp (vector times vector plus vector times vector):
Definition: Vmath.hpp:439
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.hpp:825