Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 using namespace std;
40 
41 namespace Nektar
42 {
43  string ImageWarpingSystem::className =
45  "ImageWarpingSystem",
46  ImageWarpingSystem::create,
47  "Image warping system.");
48 
49  ImageWarpingSystem::ImageWarpingSystem(
51  : UnsteadySystem(pSession)
52  {
53  }
54 
56  {
58 
59  // Define Velocity fields
61  int nq = m_fields[0]->GetNpoints();
62 
63  for(int i = 0; i < m_spacedim; ++i)
64  {
65  m_velocity[i] = Array<OneD, NekDouble> (nq,0.0);
66  }
67 
68  // Bit of a hack: redefine u/v fields so they are continuous for
69  // Helmholtz solve.
73  m_fields[2] = fld;
75  ::AllocateSharedPtr(*fld,m_graph,m_session->GetVariable(3));
76 
77  // Tell UnsteadySystem to only integrate first two fields (i.e. I and
78  // phi).
79  m_intVariables.push_back(0);
80  m_intVariables.push_back(1);
81 
83  {
86  }
87  else
88  {
89  ASSERTL0(false, "Implicit unsteady Advection not set up.");
90  }
91  }
92 
94  {
95 
96  }
97 
99  const Array<OneD, const Array<OneD,NekDouble> > &inarray,
100  Array<OneD, Array<OneD,NekDouble> > &outarray,
101  const NekDouble time)
102  {
103  int i;
104  int npoints = GetNpoints();
105  int ncoeffs = inarray[0].num_elements();
107 
108  // Load parameter alpha.
109  m_session->LoadParameter("Alpha", m_alpha);
110 
112  "CG not implemented yet.");
113 
114  // Set up storage arrays.
115  Array<OneD, NekDouble> alloc(3*npoints);
116  Array<OneD, NekDouble> dIdx1(alloc);
117  Array<OneD, NekDouble> dIdx2(alloc+npoints);
118  Array<OneD, NekDouble> dIdx3(alloc+2*npoints);
121 
122  // Calculate grad I.
123  m_fields[0]->PhysDeriv(inarray[0], dIdx1, dIdx2);
124 
125  // Set factors.
126  // TODO: Check - should be -1?
127  factors[StdRegions::eFactorLambda] = 1.0 / m_alpha / m_alpha;
128 
129  // Multiply by phi, and perform Helmholtz solve to calculate the
130  // advection velocity field.
131  for (i = 0; i < 2; ++i)
132  {
133  Vmath::Vmul(npoints, &alloc[i*npoints], 1, inarray[1].get(), 1,
134  m_fields[i+2]->UpdatePhys().get(), 1);
135  Vmath::Smul(npoints, 1/m_alpha/m_alpha, m_fields[i+2]->GetPhys().get(), 1,
136  m_fields[i+2]->UpdatePhys().get(), 1);
137  m_fields[i+2]->HelmSolve(m_fields[i+2]->GetPhys(),
138  m_fields[i+2]->UpdateCoeffs(),
139  NullFlagList, factors);
140  m_fields[i+2]->BwdTrans(m_fields[i+2]->GetCoeffs(),
141  m_velocity[i]);
142  }
143 
144  // Set up arrays for weak advection output.
145  WeakAdv[0] = Array<OneD, NekDouble>(2*ncoeffs);
146  WeakAdv[1] = WeakAdv[0]+ncoeffs;
147 
148  // Calculate the weak advection operator for I and phi - result is put
149  // in WeakAdv and is in coefficient space.
150  WeakDGAdvection(inarray, WeakAdv, true, true, 2);
151 
152  // Calculate du/dx -> dIdx1, dv/dy -> dIdx2.
153  m_fields[2]->PhysDeriv(m_velocity[0], dIdx1, dIdx3);
154  m_fields[3]->PhysDeriv(m_velocity[1], dIdx3, dIdx2);
155 
156  // Calculate RHS = I*div(u) = I*du/dx + I*dv/dy -> dIdx1.
157  Vmath::Vvtvvtp(npoints, dIdx1.get(), 1, inarray[0].get(), 1,
158  dIdx2.get(), 1, inarray[0].get(), 1,
159  dIdx1.get(), 1);
160 
161  // Take inner product to get to coefficient space.
162  Array<OneD, NekDouble> tmp2(ncoeffs);
163  m_fields[0]->IProductWRTBase(dIdx1, tmp2);
164 
165  // Add this to the weak advection for intensity field
166  // equation.
167  Vmath::Vsub(npoints, WeakAdv[0], 1, tmp2, 1, WeakAdv[0], 1);
168 
169  // Multiply by elemental inverse mass matrix, backwards transform and
170  // negate (to put on RHS of ODE).
171  for(i = 0; i < 2; ++i)
172  {
173  m_fields[i]->MultiplyByElmtInvMass(WeakAdv[i], WeakAdv[i]);
174  m_fields[i]->BwdTrans(WeakAdv[i],outarray[i]);
175  Vmath::Neg(npoints,outarray[i],1);
176  }
177  }
178 
179 
180 
181  /**
182  *
183  */
185  const Array<OneD, NekDouble> >&inarray,
186  Array<OneD, Array<OneD, NekDouble> >&outarray,
187  const NekDouble time)
188  {
189  int i;
190  int nvariables = inarray.num_elements();
191  SetBoundaryConditions(time);
192 
193  switch(m_projectionType)
194  {
196  {
197  // Just copy over array
198  int npoints = GetNpoints();
199 
200  for(i = 0; i < nvariables; ++i)
201  {
202  Vmath::Vcopy(npoints,inarray[i],1,outarray[i],1);
203  }
204  }
205  break;
206  default:
207  ASSERTL0(false,"Unknown projection scheme");
208  break;
209  }
210  }
211 
212 
214  const int i,
215  Array<OneD, Array<OneD, NekDouble> > &physfield,
217  {
218  for(int j = 0; j < flux.num_elements(); ++j)
219  {
220  Vmath::Vmul(GetNpoints(),physfield[i],1,
221  m_velocity[j],1,flux[j],1);
222  }
223  }
224 
226  Array<OneD, Array<OneD, NekDouble> > &physfield,
227  Array<OneD, Array<OneD, NekDouble> > &numflux)
228  {
229  int i;
230 
231  int nTraceNumPoints = GetTraceNpoints();
232  int nvel = m_spacedim; //m_velocity.num_elements();
233 
234  Array<OneD, NekDouble > Fwd(nTraceNumPoints);
235  Array<OneD, NekDouble > Bwd(nTraceNumPoints);
236  Array<OneD, NekDouble > Vn (nTraceNumPoints,0.0);
237 
238  // Get Edge Velocity
239  for(i = 0; i < nvel; ++i)
240  {
241  m_fields[0]->ExtractTracePhys(m_velocity[i], Fwd);
242  Vmath::Vvtvp(nTraceNumPoints,m_traceNormals[i],1,Fwd,1,Vn,1,Vn,1);
243  }
244 
245  for(i = 0; i < numflux.num_elements(); ++i)
246  {
247  m_fields[i]->GetFwdBwdTracePhys(physfield[i],Fwd,Bwd);
248  // Evaulate upwinded m_fields[i]
249  m_fields[i]->GetTrace()->Upwind(Vn,Fwd,Bwd,numflux[i]);
250  // Calculate m_fields[i]*Vn
251  Vmath::Vmul(nTraceNumPoints,numflux[i],1,Vn,1,numflux[i],1);
252  }
253  }
254 
256  {
258  }
259 }
Array< OneD, Array< OneD, NekDouble > > m_velocity
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
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:442
enum MultiRegions::ProjectionType m_projectionType
Type of projection; e.g continuous or discontinuous.
STL namespace.
boost::shared_ptr< ContField2D > ContField2DSharedPtr
Definition: ContField2D.h:289
std::map< ConstFactorType, NekDouble > ConstFactorMap
Definition: StdRegions.hpp:252
boost::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: MeshPartition.h:51
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)
virtual SOLVER_UTILS_EXPORT void v_GenerateSummary(SummaryList &s)
Print a summary of time stepping parameters.
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:213
bool m_explicitAdvection
Indicates if explicit or implicit treatment of advection is used.
void DefineOdeRhs(FuncPointerT func, ObjectPointerT obj)
Base class for unsteady solvers.
int m_spacedim
Spatial dimension (>= expansion dim).
void Neg(int n, T *x, const int incx)
Negate x = -x.
Definition: Vmath.cpp:396
double NekDouble
virtual SOLVER_UTILS_EXPORT void v_InitObject()
Init object for UnsteadySystem class.
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:343
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:537
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.
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:1061
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:183
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