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 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Image warping solve routines
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
38 
39 namespace Nektar
40 {
43  "ImageWarpingSystem",
45  "Image warping system.");
46 
49  : UnsteadySystem(pSession)
50  {
51  }
52 
54  {
55  UnsteadySystem::v_InitObject();
56 
57  // Define Velocity fields
59  int nq = m_fields[0]->GetNpoints();
60 
61  for(int i = 0; i < m_spacedim; ++i)
62  {
63  m_velocity[i] = Array<OneD, NekDouble> (nq,0.0);
64  }
65 
66  // Bit of a hack: redefine u/v fields so they are continuous for
67  // Helmholtz solve.
71  m_fields[2] = fld;
73  ::AllocateSharedPtr(*fld,m_graph,m_session->GetVariable(3));
74 
75  // Tell UnsteadySystem to only integrate first two fields (i.e. I and
76  // phi).
77  m_intVariables.push_back(0);
78  m_intVariables.push_back(1);
79 
81  {
84  }
85  else
86  {
87  ASSERTL0(false, "Implicit unsteady Advection not set up.");
88  }
89  }
90 
92  {
93 
94  }
95 
97  const Array<OneD, const Array<OneD,NekDouble> > &inarray,
98  Array<OneD, Array<OneD,NekDouble> > &outarray,
99  const NekDouble time)
100  {
101  int i;
102  int npoints = GetNpoints();
103  int ncoeffs = inarray[0].num_elements();
105 
106  // Load parameter alpha.
107  m_session->LoadParameter("Alpha", m_alpha);
108 
110  "CG not implemented yet.");
111 
112  // Set up storage arrays.
113  Array<OneD, NekDouble> alloc(3*npoints);
114  Array<OneD, NekDouble> dIdx1(alloc);
115  Array<OneD, NekDouble> dIdx2(alloc+npoints);
116  Array<OneD, NekDouble> dIdx3(alloc+2*npoints);
119 
120  // Calculate grad I.
121  m_fields[0]->PhysDeriv(inarray[0], dIdx1, dIdx2);
122 
123  // Set factors.
124  // TODO: Check - should be -1?
125  factors[StdRegions::eFactorLambda] = 1.0 / m_alpha / m_alpha;
126 
127  // Multiply by phi, and perform Helmholtz solve to calculate the
128  // advection velocity field.
129  for (i = 0; i < 2; ++i)
130  {
131  Vmath::Vmul(npoints, &alloc[i*npoints], 1, inarray[1].get(), 1,
132  m_fields[i+2]->UpdatePhys().get(), 1);
133  Vmath::Smul(npoints, 1/m_alpha/m_alpha, m_fields[i+2]->GetPhys().get(), 1,
134  m_fields[i+2]->UpdatePhys().get(), 1);
135  m_fields[i+2]->HelmSolve(m_fields[i+2]->GetPhys(),
136  m_fields[i+2]->UpdateCoeffs(),
137  NullFlagList, factors);
138  m_fields[i+2]->BwdTrans(m_fields[i+2]->GetCoeffs(),
139  m_velocity[i]);
140  }
141 
142  // Set up arrays for weak advection output.
143  WeakAdv[0] = Array<OneD, NekDouble>(2*ncoeffs);
144  WeakAdv[1] = WeakAdv[0]+ncoeffs;
145 
146  // Calculate the weak advection operator for I and phi - result is put
147  // in WeakAdv and is in coefficient space.
148  WeakDGAdvection(inarray, WeakAdv, true, true, 2);
149 
150  // Calculate du/dx -> dIdx1, dv/dy -> dIdx2.
151  m_fields[2]->PhysDeriv(m_velocity[0], dIdx1, dIdx3);
152  m_fields[3]->PhysDeriv(m_velocity[1], dIdx3, dIdx2);
153 
154  // Calculate RHS = I*div(u) = I*du/dx + I*dv/dy -> dIdx1.
155  Vmath::Vvtvvtp(npoints, dIdx1.get(), 1, inarray[0].get(), 1,
156  dIdx2.get(), 1, inarray[0].get(), 1,
157  dIdx1.get(), 1);
158 
159  // Take inner product to get to coefficient space.
160  Array<OneD, NekDouble> tmp2(ncoeffs);
161  m_fields[0]->IProductWRTBase(dIdx1, tmp2);
162 
163  // Add this to the weak advection for intensity field
164  // equation.
165  Vmath::Vsub(npoints, WeakAdv[0], 1, tmp2, 1, WeakAdv[0], 1);
166 
167  // Multiply by elemental inverse mass matrix, backwards transform and
168  // negate (to put on RHS of ODE).
169  for(i = 0; i < 2; ++i)
170  {
171  m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
172  m_fields[i]->BwdTrans(WeakAdv[i],outarray[i]);
173  Vmath::Neg(npoints,outarray[i],1);
174  }
175  }
176 
177 
178 
179  /**
180  *
181  */
183  const Array<OneD, NekDouble> >&inarray,
184  Array<OneD, Array<OneD, NekDouble> >&outarray,
185  const NekDouble time)
186  {
187  int i;
188  int nvariables = inarray.num_elements();
189  SetBoundaryConditions(time);
190 
191  switch(m_projectionType)
192  {
194  {
195  // Just copy over array
196  int npoints = GetNpoints();
197 
198  for(i = 0; i < nvariables; ++i)
199  {
200  Vmath::Vcopy(npoints,inarray[i],1,outarray[i],1);
201  }
202  }
203  break;
204  default:
205  ASSERTL0(false,"Unknown projection scheme");
206  break;
207  }
208  }
209 
210 
212  const int i,
213  Array<OneD, Array<OneD, NekDouble> > &physfield,
215  {
216  for(int j = 0; j < flux.num_elements(); ++j)
217  {
218  Vmath::Vmul(GetNpoints(),physfield[i],1,
219  m_velocity[j],1,flux[j],1);
220  }
221  }
222 
224  Array<OneD, Array<OneD, NekDouble> > &physfield,
225  Array<OneD, Array<OneD, NekDouble> > &numflux)
226  {
227  int i;
228 
229  int nTraceNumPoints = GetTraceNpoints();
230  int nvel = m_spacedim; //m_velocity.num_elements();
231 
232  Array<OneD, NekDouble > Fwd(nTraceNumPoints);
233  Array<OneD, NekDouble > Bwd(nTraceNumPoints);
234  Array<OneD, NekDouble > Vn (nTraceNumPoints,0.0);
235 
236  // Get Edge Velocity
237  for(i = 0; i < nvel; ++i)
238  {
239  m_fields[0]->ExtractTracePhys(m_velocity[i], Fwd);
240  Vmath::Vvtvp(nTraceNumPoints,m_traceNormals[i],1,Fwd,1,Vn,1,Vn,1);
241  }
242 
243  for(i = 0; i < numflux.num_elements(); ++i)
244  {
245  m_fields[i]->GetFwdBwdTracePhys(physfield[i],Fwd,Bwd);
246  // Evaulate upwinded m_fields[i]
247  m_fields[i]->GetTrace()->Upwind(Vn,Fwd,Bwd,numflux[i]);
248  // Calculate m_fields[i]*Vn
249  Vmath::Vmul(nTraceNumPoints,numflux[i],1,Vn,1,numflux[i],1);
250  }
251  }
252 
254  {
255  UnsteadySystem::v_GenerateSummary(s);
256  }
257 }
Array< OneD, Array< OneD, NekDouble > > m_velocity
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:161
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
void DoOdeRhs(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
LibUtilities::TimeIntegrationSchemeOperators m_ode
The time integration scheme operators to use.
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:47
virtual void v_NumericalFlux(Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &numflux)
virtual void v_InitObject()
Init object for UnsteadySystem class.
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.cpp:428
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
static EquationSystemSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession)
Creates an instance of this class.
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:291
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:248
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:50
virtual void v_GetFluxVector(const int i, Array< OneD, Array< OneD, NekDouble > > &physfield, Array< OneD, Array< OneD, NekDouble > > &flux)
Array< OneD, Array< OneD, NekDouble > > m_traceNormals
Array holding trace normals for DG simulations in the forwards direction.
virtual void v_GenerateSummary(SolverUtils::SummaryList &s)
Print a summary of time stepping parameters.
void DefineProjection(FuncPointerT func, ObjectPointerT obj)
void Smul(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Scalar multiply y = alpha*y.
Definition: Vmath.cpp:199
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
static std::string className
Name of class.
int m_spacedim
Spatial dimension (>= expansion dim).
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:382
double NekDouble
EquationSystemFactory & GetEquationSystemFactory()
void DoOdeProjection(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
SOLVER_UTILS_EXPORT void SetBoundaryConditions(NekDouble time)
Evaluates the boundary conditions at the given time.
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.cpp:329
SOLVER_UTILS_EXPORT int GetNpoints()
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.cpp:523
Array< OneD, MultiRegions::ExpListSharedPtr > m_fields
Array holding all dependent variables.
LibUtilities::SessionReaderSharedPtr m_session
The session reader.
SOLVER_UTILS_EXPORT int GetTraceNpoints()
SpatialDomains::MeshGraphSharedPtr m_graph
Pointer to graph defining mesh.
ImageWarpingSystem(const LibUtilities::SessionReaderSharedPtr &pSession)
SOLVER_UTILS_EXPORT void WeakDGAdvection(const Array< OneD, Array< OneD, NekDouble > > &InField, Array< OneD, Array< OneD, NekDouble > > &OutField, bool NumericalFluxIncludesNormal=true, bool InFieldIsInPhysSpace=false, int nvariables=0)
Calculate the weak discontinuous Galerkin advection.
void Vcopy(int n, const T *x, const int incx, T *y, const int incy)
Definition: Vmath.cpp:1038
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.cpp:169
static FlagList NullFlagList
An empty flag list.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215