Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ProcessCurve.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: ProcessCurve.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: create curved edges using a custom expression y = f(x)
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #include <LocalRegions/SegExp.h>
37 #include <LocalRegions/QuadExp.h>
38 #include <LocalRegions/TriExp.h>
40 
43 
45 
47 
48 #include "ProcessCurve.h"
49 
50 using namespace std;
51 using namespace Nektar::NekMeshUtils;
52 
53 namespace Nektar
54 {
55 namespace Utilities
56 {
57 
58 ModuleKey ProcessCurve::className = GetModuleFactory().RegisterCreatorFunction(
59  ModuleKey(eProcessModule, "curve"),
60  ProcessCurve::create,
61  "Creates curved edge information using a function y=f(x) (2D only).");
62 
63 /**
64  * @brief Default constructor.
65  */
66 ProcessCurve::ProcessCurve(MeshSharedPtr m) : ProcessCurvedEdges(m)
67 {
68  m_config["function"] = ConfigOption(false, "NotSet",
69  "Expression of the curve: y = f(x).");
70  m_config["file"] = ConfigOption(false, "NotSet",
71  "Pts file containing coordinates (x,y).");
72  m_config["niter"] = ConfigOption(false, "50",
73  "Number of iterations to perform to obtain evenly distribution of points.");
74  m_config["gamma"] = ConfigOption(false, "0.1",
75  "Relaxation parameter.");
76 }
77 
78 /**
79  * @brief Destructor.
80  */
82 {
83 }
84 
85 /**
86  * This function generates curved edge information in 2D, by including
87  * equispaced nodes to an edge following the prescribed function \f$y = f(x)\f$.
88  *
89  * An approximately equispaced node distribution is obtained by initially
90  * considering a uniform distribution in the \f$ x\f$ direction, and then
91  * adjusting the interior nodes using the iteration
92  * \f[ x^{n}_{i+1} = x^{n-1}_{i+1} + \Delta x^{n}_{i+1} \f]
93  * \f[ \Delta x^{n}_{i+1} = \Delta x^{n}_{i} +
94  * \gamma (\overline{\Delta s} - \Delta s^{n}_{i}) sign(\Delta x^{n}_{i}) \f]
95  *
96  * where \f$ x^{n}\f$ is the \f$ x\f$ coordinate of the \f$ n^{th}\f$ node,
97  * \f$ i\f$ is the iteration counter, \f$ \Delta s^{n}\f$ is the Cartesian
98  * distance between nodes \f$ n\f$ and \f$ n-1\f$ and \f$ \overline{\Delta s}\f$
99  * is the average of these distances. The relaxation factor \f$ \gamma\f$ and
100  * the fixed number of iterations \f$ Niter\f$ are parameters of the module.
101  *
102  * In case the correction to \f$ \Delta x\f$ leads to an invalid distribution
103  * (i.e. reversing directions), the relaxation parameter is successively halved
104  * until it leads to a valid distribution.
105  *
106  * @param edge Edge which will be modified
107  */
109 {
110  NodeSharedPtr n1 = edge->m_n1;
111  NodeSharedPtr n2 = edge->m_n2;
112 
113  int nq = m_config["N"].as<int>();
114  int niter = m_config["niter"].as<int>();
115  NekDouble gamma = m_config["gamma"].as<double>();
116 
117  edge->m_edgeNodes.resize(nq - 2);
118 
119  // Read function defining the curve
120  if (m_config["function"].as<string>().compare("NotSet") != 0)
121  {
122  ASSERTL0(m_config["file"].as<string>().compare("NotSet") == 0,
123  "Function and file cannot be defined at the same time.");
124 
125  std::string fstr = m_config["function"].as<string>();
126  m_fExprId = m_fEval.DefineFunction("x y z", fstr);
127  m_fromFile = false;
128  }
129  else
130  {
131  ASSERTL0(m_config["file"].as<string>().compare("NotSet") != 0,
132  "Need to define either function or file.");
133  std::string inFile = m_config["file"].as<string>();
134 
139  ptsIO->Import(inFile, m_fieldPts);
140 
141  m_fromFile = true;
142  }
143 
144  // Coordinates of points
145  Array<OneD, NekDouble> x(nq,0.0);
146  Array<OneD, NekDouble> y(nq,0.0);
147  // Distances between points
148  Array<OneD, NekDouble> dx(nq-1,0.0);
149  Array<OneD, NekDouble> dy(nq-1,0.0);
150  Array<OneD, NekDouble> ds(nq-1,0.0);
151  // Average distance and deviation from average
152  Array<OneD, NekDouble> s_deviation(nq-1,0.0);
153  NekDouble s_average;
154 
155  // Fix start point
156  x[0] = n1->m_x;
157  y[0] = EvaluateCoordinate(x[0]);
158  // Start with uniform distribution along x-axis
159  Vmath::Sadd(nq-1, (n2->m_x - n1->m_x) / (nq-1), dx, 1, dx, 1);
160 
161  // Iterate a few times to make points more evenly distributed
162  for (int s = 0; s < niter; ++s)
163  {
164  s_average = 0.0;
165  for (int k = 1; k < nq; ++k)
166  {
167  x[k] = x[k-1] + dx[k-1] + gamma*s_deviation[k-1];
168  y[k] = EvaluateCoordinate(x[k]);
169 
170  dx[k-1] = x[k] - x[k-1];
171  dy[k-1] = y[k] - y[k-1];
172  ds[k-1] = sqrt(dx[k-1]*dx[k-1] + dy[k-1]*dy[k-1]);
173 
174  s_average = s_average + ds[k-1]/(nq-1);
175  }
176  // Calculate correction for next iteration
177  for (int k = 0; k < nq-1; ++k)
178  {
179  s_deviation[k] = (s_average - ds[k])*
180  ( dx[k]/abs(dx[k]));
181  }
182  // Adjust gama to make sure next partition is valid
183  // (no reversals in dx)
184  bool valid = false;
185  gamma = m_config["gamma"].as<double>();
186  while( !valid)
187  {
188  valid = true;
189  for (int k = 0; k < nq-1; ++k)
190  {
191  if( dx[k]*(dx[k]+gamma*s_deviation[k]) < 0.0)
192  {
193  gamma = gamma/2;
194  valid = false;
195  continue;
196  }
197  }
198  }
199  }
200 
201  // Write interior nodes to edge
202  for (int k = 1; k < nq-1; ++k)
203  {
204  edge->m_edgeNodes[k-1] = NodeSharedPtr(
205  new Node(0, x[k], y[k], n1->m_z));
206  }
207  edge->m_curveType = LibUtilities::ePolyEvenlySpaced;
208 }
209 
211 {
212  if (m_fromFile)
213  {
215  tmp[0] = Array<OneD, NekDouble>(1, xCoord);
216  tmp[1] = Array<OneD, NekDouble>(1, 0.0);
217 
220 
222  interp.Interpolate(m_fieldPts, toPts);
223 
224  return tmp[1][0];
225  }
226  else
227  {
228  return m_fEval.Evaluate(m_fExprId, xCoord, 0.0, 0.0, 0.0);
229  }
230 }
231 
232 }
233 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
A class that contains algorithms for interpolation between pts fields, expansions and different meshe...
Definition: Interpolator.h:78
tBaseSharedPtr CreateInstance(tKey idKey BOOST_PP_COMMA_IF(MAX_PARAM) BOOST_PP_ENUM_BINARY_PARAMS(MAX_PARAM, tParam, x))
Create an instance of the class referred to by idKey.
Definition: NekFactory.hpp:162
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
FIELD_UTILS_EXPORT void Interpolate(const LibUtilities::PtsFieldSharedPtr ptsInField, LibUtilities::PtsFieldSharedPtr &ptsOutField)
Interpolate from a pts field to a pts field.
LibUtilities::AnalyticExpressionEvaluator m_fEval
Definition: ProcessCurve.h:68
NekDouble Evaluate(const int AnalyticExpression_id)
Evaluation method for expressions depending on parameters only.
STL namespace.
pair< ModuleType, string > ModuleKey
CommFactory & GetCommFactory()
Definition: Comm.cpp:61
boost::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:178
1D Evenly-spaced points using Lagrange polynomial
Definition: PointsType.h:65
virtual ~ProcessCurve()
Destructor.
LibUtilities::PtsFieldSharedPtr m_fieldPts
Definition: ProcessCurve.h:70
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
Represents a command-line configuration option.
boost::shared_ptr< Node > NodeSharedPtr
Definition: Node.h:50
double NekDouble
std::map< std::string, ConfigOption > m_config
List of configuration values.
void Sadd(int n, const T alpha, const T *x, const int incx, T *y, const int incy)
Add vector y = alpha + x.
Definition: Vmath.cpp:315
boost::shared_ptr< PtsIO > PtsIOSharedPtr
Definition: PtsIO.h:95
void v_GenerateEdgeNodes(NekMeshUtils::EdgeSharedPtr edge)
boost::shared_ptr< Edge > EdgeSharedPtr
Shared pointer to an edge.
Definition: Edge.h:135
int DefineFunction(const std::string &vlist, const std::string &function)
This function allows one to define a function to evaluate. The first argument (vlist) is a list of va...
boost::shared_ptr< Mesh > MeshSharedPtr
Shared pointer to a mesh.
Definition: Mesh.h:147
NekDouble EvaluateCoordinate(NekDouble xCoord)
std::pair< ModuleType, std::string > ModuleKey
ModuleFactory & GetModuleFactory()
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215