young_laplace_elements.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 
27 // Non-inline functions for YoungLaplace elements
28 #include "young_laplace_elements.h"
29 
30 namespace oomph
31 {
32  //======================================================================
33  // Set the data for the number of Variables at each node
34  //======================================================================
35  template<>
37  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
38  template<>
40  1, 1, 1, 1, 1, 1, 1, 1, 1};
41  template<>
42  const unsigned QYoungLaplaceElement<2>::Initial_Nvalue[4] = {1, 1, 1, 1};
43 
44 
45  //======================================================================
46  /// Get exact position vector to meniscus
47  //======================================================================
49  const Vector<double>& s,
50  Vector<double>& r,
52  {
53  // Get global coordinates
54  Vector<double> x(2);
55  interpolated_x(s, x);
56 
57  // Exact solution Vector (here a scalar)
58  Vector<double> exact_soln(1);
59 
60  // Get exact solution at this point
61  (*exact_soln_pt)(x, exact_soln);
62 
63  if (!use_spines())
64  {
65  r[0] = x[0];
66  r[1] = x[1];
67  r[2] = exact_soln[0];
68  }
69  else
70  {
71  /// Get spines values
72  Vector<double> spine_base(3, 0.0);
73  Vector<double> spine(3, 0.0);
74  Vector<Vector<double>> dspine_base;
75  allocate_vector_of_vectors(2, 3, dspine_base);
76  Vector<Vector<double>> dspine;
77  allocate_vector_of_vectors(2, 3, dspine);
78 
79  get_spine_base(x, spine_base, dspine_base);
80  get_spine(x, spine, dspine);
81 
82  /// Global Eulerian cooordinates
83  for (unsigned j = 0; j < 3; j++)
84  {
85  r[j] = spine_base[j] + exact_soln[0] * spine[j];
86  }
87  }
88  }
89 
90  //======================================================================
91  /// Get position vector to meniscus
92  //======================================================================
94  Vector<double>& r) const
95  {
96  // Get global coordinates
97  Vector<double> x(2);
98  interpolated_x(s, x);
99 
100  // Displacement along spine (or cartesian displacement)
101  double u = interpolated_u(s);
102 
103  // cartesian calculation case
104  if (!use_spines())
105  {
106  r[0] = x[0];
107  r[1] = x[1];
108  r[2] = u;
109  }
110  // spine case
111  else
112  {
113  /// Get spines values
114  Vector<double> spine_base(3, 0.0);
115  Vector<double> spine(3, 0.0);
116  Vector<Vector<double>> dspine_base;
117  allocate_vector_of_vectors(2, 3, dspine_base);
118  Vector<Vector<double>> dspine;
119  allocate_vector_of_vectors(2, 3, dspine);
120 
121  get_spine_base(x, spine_base, dspine_base);
122  get_spine(x, spine, dspine);
123 
124  /// Global Eulerian cooordinates
125  for (unsigned j = 0; j < 3; j++)
126  {
127  r[j] = spine_base[j] + u * spine[j];
128  }
129  }
130  }
131 
132 
133  //======================================================================
134  /// Compute element residual vector. Pure version without hanging nodes
135  //======================================================================
137  Vector<double>& residuals)
138  {
139  // Find out how many nodes there are
140  unsigned n_node = nnode();
141 
142  // Set up memory for the shape functions
143  Shape psi(n_node);
144  DShape dpsidzeta(n_node, 2);
145 
146  // Set the value of n_intpt
147  unsigned n_intpt = integral_pt()->nweight();
148 
149  // Integers to store the local equation numbers
150  int local_eqn = 0;
151 
152  // Loop over the integration points
153  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
154  {
155  // Get the integral weight
156  double w = integral_pt()->weight(ipt);
157 
158  // Call the derivatives of the shape and test functions
159  double J = dshape_eulerian_at_knot(ipt, psi, dpsidzeta);
160 
161  // Premultiply the weights and the Jacobian
162  double W = w * J;
163 
164  // Calculate local values of displacement along spine and its derivatives
165  // Allocate and initialise to zero
166  double interpolated_u = 0.0;
168  Vector<double> interpolated_dudzeta(2, 0.0);
169 
170  // Calculate function value and derivatives:
171  //-----------------------------------------
172  // Loop over nodes
173  for (unsigned l = 0; l < n_node; l++)
174  {
175  interpolated_u += u(l) * psi(l);
176  // Loop over directions
177  for (unsigned j = 0; j < 2; j++)
178  {
179  interpolated_zeta[j] += nodal_position(l, j) * psi(l);
180  interpolated_dudzeta[j] += u(l) * dpsidzeta(l, j);
181  }
182  }
183 
184 
185  // Allocation and definition of variables necessary for
186  // further calculations
187 
188  /// "Simple" case
189  /// --------------
190  double nonlinearterm = 1.0;
191  double sqnorm = 0.0;
192 
193  /// Spine case
194  /// -----------
195 
196  // Derivs of position vector w.r.t. global intrinsic coords
197  Vector<Vector<double>> dRdzeta;
198  allocate_vector_of_vectors(2, 3, dRdzeta);
199 
200  // Unnormalised normal
201  Vector<double> N_unnormalised(3, 0.0);
202 
203  // Spine and spine basis vectors, entries initialised to zero
204  Vector<double> spine_base(3, 0.0), spine(3, 0.0);
205 
206  // Derivative of spine basis vector w.r.t to the intrinsic
207  // coordinates: dspine_base[i,j] = j-th component of the deriv.
208  // of the spine basis vector w.r.t. to the i-th global intrinsic
209  // coordinate
210  Vector<Vector<double>> dspine_base;
211  allocate_vector_of_vectors(2, 3, dspine_base);
212 
213  // Derivative of spine vector w.r.t to the intrinsic
214  // coordinates: dspine[i,j] = j-th component of the deriv.
215  // of the spine vector w.r.t. to the i-th global intrinsic
216  // coordinate
217  Vector<Vector<double>> dspine;
218  allocate_vector_of_vectors(2, 3, dspine);
219 
220  // Vector v_\alpha contains the numerator of the variations of the
221  // area element {\cal A}^{1/2} w.r.t. the components of dR/d\zeta_\alpha
222  Vector<double> area_variation_numerator_0(3, 0.0);
223  Vector<double> area_variation_numerator_1(3, 0.0);
224 
225  // Vector position
226  Vector<double> r(3, 0.0);
227 
228  // No spines
229  //---------
230  if (!use_spines())
231  {
232  for (unsigned j = 0; j < 2; j++)
233  {
234  sqnorm += interpolated_dudzeta[j] * interpolated_dudzeta[j];
235  }
236  nonlinearterm = 1.0 / sqrt(1.0 + sqnorm);
237  }
238 
239  // Spines
240  //------
241  else
242  {
243  // Get the spines
244  get_spine_base(interpolated_zeta, spine_base, dspine_base);
245  get_spine(interpolated_zeta, spine, dspine);
246 
247  // calculation of dR/d\zeta_\alpha
248  for (unsigned alpha = 0; alpha < 2; alpha++)
249  {
250  // Product rule for d(u {\bf S} ) / d \zeta_\alpha
251  Vector<double> dudzeta_times_spine(3, 0.0);
253  interpolated_dudzeta[alpha], spine, dudzeta_times_spine);
254 
255  Vector<double> u_times_dspinedzeta(3, 0.0);
257  interpolated_u, dspine[alpha], u_times_dspinedzeta);
258 
259  Vector<double> d_u_times_spine_dzeta(3, 0.0);
260  vector_sum(
261  dudzeta_times_spine, u_times_dspinedzeta, d_u_times_spine_dzeta);
262 
263  // Add derivative of spine base
264  vector_sum(d_u_times_spine_dzeta, dspine_base[alpha], dRdzeta[alpha]);
265  }
266 
267  /// Get the unnormalized normal
268  cross_product(dRdzeta[0], dRdzeta[1], N_unnormalised);
269 
270  // Tmp storage
271  Vector<double> v_tmp_1(3, 0.0);
272  Vector<double> v_tmp_2(3, 0.0);
273 
274  // Calculation of
275  // |dR/d\zeta_1|^2 dR/d\zeta_0 - <dR/d\zeta_0,dR/d\zeta_1>dR/d\zeta_1
276  scalar_times_vector(pow(two_norm(dRdzeta[1]), 2), dRdzeta[0], v_tmp_1);
278  -1 * scalar_product(dRdzeta[0], dRdzeta[1]), dRdzeta[1], v_tmp_2);
279  vector_sum(v_tmp_1, v_tmp_2, area_variation_numerator_0);
280 
281  // Calculation of
282  // |dR/d\zeta_0|^2 dR/d\zeta_1 - <dR/d\zeta_0,dR/d\zeta_1>dR/d\zeta_0
283  scalar_times_vector(pow(two_norm(dRdzeta[0]), 2), dRdzeta[1], v_tmp_1);
285  -1 * scalar_product(dRdzeta[0], dRdzeta[1]), dRdzeta[0], v_tmp_2);
286  vector_sum(v_tmp_1, v_tmp_2, area_variation_numerator_1);
287 
288  // Global Eulerian cooordinates
289  for (unsigned j = 0; j < 3; j++)
290  {
291  r[j] = spine_base[j] + interpolated_u * spine[j];
292  }
293  }
294 
295 
296  // Assemble residuals
297  //-------------------
298 
299  // Loop over the test (shape) functions
300  for (unsigned l = 0; l < n_node; l++)
301  {
302  // Get the local equation
303  local_eqn = u_local_eqn(l);
304 
305  /*IF it's not a boundary condition*/
306  if (local_eqn >= 0)
307  {
308  // "simple" calculation case
309  if (!use_spines())
310  {
311  // Add source term: The curvature
312  residuals[local_eqn] += get_kappa() * psi(l) * W;
313 
314  // The YoungLaplace bit itself
315  for (unsigned k = 0; k < 2; k++)
316  {
317  residuals[local_eqn] +=
318  nonlinearterm * interpolated_dudzeta[k] * dpsidzeta(l, k) * W;
319  }
320  }
321 
322  // Spine calculation case
323  else
324  {
325  // Calculation of d(u S)/d\zeta_0
326  //-------------------------------
327  Vector<double> v_tmp_1(3, 0.0);
328  scalar_times_vector(dpsidzeta(l, 0), spine, v_tmp_1);
329 
330  Vector<double> v_tmp_2(3, 0.0);
331  scalar_times_vector(psi(l), dspine[0], v_tmp_2);
332 
333  Vector<double> d_uS_dzeta0(3, 0.0);
334  vector_sum(v_tmp_1, v_tmp_2, d_uS_dzeta0);
335 
336  // Add contribution to residual
337  residuals[local_eqn] +=
338  W * scalar_product(area_variation_numerator_0, d_uS_dzeta0) /
339  two_norm(N_unnormalised);
340 
341  // Calculation of d(u S)/d\zeta_1
342  scalar_times_vector(dpsidzeta(l, 1), spine, v_tmp_1);
343  scalar_times_vector(psi(l), dspine[1], v_tmp_2);
344  Vector<double> d_uS_dzeta1(3, 0.0);
345  vector_sum(v_tmp_1, v_tmp_2, d_uS_dzeta1);
346 
347  // Add contribution to residual
348  residuals[local_eqn] +=
349  W * scalar_product(area_variation_numerator_1, d_uS_dzeta1) /
350  two_norm(N_unnormalised);
351 
352  // Curvature contribution to the residual : kappa N S test
353  residuals[local_eqn] += W * (get_kappa()) *
354  scalar_product(N_unnormalised, spine) *
355  psi(l);
356  }
357  }
358  }
359 
360  } // End of loop over integration points
361  }
362 
363 
364  //======================================================================
365  /// Self-test: Return 0 for OK
366  //======================================================================
368  {
369  bool passed = true;
370 
371  // Check lower-level stuff
372  if (FiniteElement::self_test() != 0)
373  {
374  passed = false;
375  }
376 
377  // Return verdict
378  if (passed)
379  {
380  return 0;
381  }
382  else
383  {
384  return 1;
385  }
386  }
387 
388 
389  //======================================================================
390  /// Output solution at nplot points in each coordinate direction
391  //======================================================================
392  void YoungLaplaceEquations::output(std::ostream& outfile,
393  const unsigned& nplot)
394  {
395  // Vector of local coordinates
396  Vector<double> s(2);
397 
398  // Tecplot header info
399  outfile << tecplot_zone_string(nplot);
400 
401  // Loop over plot points
402  unsigned num_plot_points = nplot_points(nplot);
403  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
404  {
405  // Get local coordinates of plot point
406  get_s_plot(iplot, nplot, s);
407 
408  // Compute intrinsic coordinates
409  Vector<double> xx(2, 0.0);
410  for (unsigned i = 0; i < 2; i++)
411  {
412  xx[i] = interpolated_x(s, i);
413  }
414 
415  // Calculate the cartesian coordinates of point on meniscus
416  Vector<double> r(3, 0.0);
417 
418  // Position
419  if (use_spines())
420  {
421  position(s, r);
422  }
423  else
424  {
425  r[0] = xx[0];
426  r[1] = xx[1];
427  r[2] = interpolated_u(s);
428  }
429 
430  // Output positon on meniscus
431  for (unsigned i = 0; i < 3; i++)
432  {
433  outfile << r[i] << " ";
434  }
435 
436  // Get spine stuff
437  Vector<double> spine_base(3, 0.0), spine(3, 0.0);
438  Vector<Vector<double>> dspine_base;
439  allocate_vector_of_vectors(2, 3, dspine_base);
440  Vector<Vector<double>> dspine;
441  allocate_vector_of_vectors(2, 3, dspine);
442 
443  // Get the spines
444  if (use_spines())
445  {
446  get_spine_base(xx, spine_base, dspine_base);
447  get_spine(xx, spine, dspine);
448  }
449 
450 
451  // Output spine base
452  for (unsigned i = 0; i < 3; i++)
453  {
454  outfile << spine_base[i] << " ";
455  }
456 
457  // Output spines
458  for (unsigned i = 0; i < 3; i++)
459  {
460  outfile << spine[i] << " ";
461  }
462 
463 
464  // Output intrinsic coordinates
465  for (unsigned i = 0; i < 2; i++)
466  {
467  outfile << xx[i] << " ";
468  }
469 
470  // Output unknown
471  outfile << interpolated_u(s) << " ";
472 
473 
474  // Done
475  outfile << std::endl;
476  }
477 
478  // Write tecplot footer (e.g. FE connectivity lists)
479  write_tecplot_zone_footer(outfile, nplot);
480  }
481 
482 
483  //======================================================================
484  /// Output exact solution
485  ///
486  /// Solution is provided via function pointer.
487  /// Plot at a given number of plot points.
488  //======================================================================
490  std::ostream& outfile,
491  const unsigned& nplot,
493  {
494  // Vector of local coordinates
495  Vector<double> s(2);
496 
497  // Vector for coordinates
498  Vector<double> x(2);
499 
500  // Tecplot header info
501  outfile << tecplot_zone_string(nplot);
502 
503  // Exact solution Vector (here a scalar)
504  Vector<double> exact_soln(1);
505 
506  // Loop over plot points
507  unsigned num_plot_points = nplot_points(nplot);
508  for (unsigned iplot = 0; iplot < num_plot_points; iplot++)
509  {
510  // Get local coordinates of plot point
511  get_s_plot(iplot, nplot, s);
512 
513  // Get x position as Vector
514  interpolated_x(s, x);
515 
516  /// Calculate the cartesian coordinates of point on meniscus
517  Vector<double> r_exact(3, 0.0);
518  exact_position(s, r_exact, exact_soln_pt);
519 
520  // Output x_exact,y_exact,z_exact
521  for (unsigned i = 0; i < 3; i++)
522  {
523  outfile << r_exact[i] << " ";
524  }
525 
526  // Done
527  outfile << std::endl;
528  }
529 
530  // Write tecplot footer (e.g. FE connectivity lists)
531  write_tecplot_zone_footer(outfile, nplot);
532  }
533 
534 
535  //======================================================================
536  /// Validate against exact solution
537  ///
538  /// Solution is provided via function pointer.
539  /// Plot error at a given number of plot points.
540  ///
541  //======================================================================
543  std::ostream& outfile,
545  double& error,
546  double& norm)
547  {
548  // Initialise
549  error = 0.0;
550  norm = 0.0;
551 
552  // Vector of local coordinates
553  Vector<double> s(2);
554 
555  // Vector for coordinates
556  Vector<double> x(2);
557 
558  // Set the value of n_intpt
559  unsigned n_intpt = integral_pt()->nweight();
560 
561  // Tecplot
562  outfile << "ZONE" << std::endl;
563 
564  // Exact solution Vector (here a scalar)
565  Vector<double> exact_soln(1);
566 
567  // Loop over the integration points
568  for (unsigned ipt = 0; ipt < n_intpt; ipt++)
569  {
570  // Assign values of s
571  for (unsigned i = 0; i < 2; i++)
572  {
573  s[i] = integral_pt()->knot(ipt, i);
574  }
575 
576  // Get the integral weight
577  double w = integral_pt()->weight(ipt);
578 
579  // Get jacobian of mapping
580  double J = J_eulerian(s);
581 
582  // Premultiply the weights and the Jacobian
583  double W = w * J;
584 
585  /// Calculate the cartesian coordinates of point on meniscus
586  Vector<double> r(3, 0.0);
587  position(s, r);
588 
589  /// Calculate the exact position
590  Vector<double> r_exact(3, 0.0);
591  exact_position(s, r_exact, exact_soln_pt);
592 
593  // Output x,y,...,error
594  for (unsigned i = 0; i < 2; i++)
595  {
596  outfile << r[i] << " ";
597  }
598 
599  for (unsigned i = 0; i < 2; i++)
600  {
601  outfile << r_exact[i] << " ";
602  }
603 
604  outfile << std::endl;
605 
606  // Add to error and norm
607  norm += 0.0;
608  for (unsigned i = 0; i < 2; i++)
609  {
610  error += (r[i] - r_exact[i]) * (r[i] - r_exact[i]) * W;
611  }
612  }
613  }
614 
615 
616  //====================================================================
617  // Force build of templates
618  //====================================================================
619  template class QYoungLaplaceElement<2>;
620  template class QYoungLaplaceElement<3>;
621  template class QYoungLaplaceElement<4>;
622 
623 } // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
Definition: elements.cc:4103
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3161
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3325
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4675
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
Definition: elements.h:3148
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
Definition: elements.h:3186
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
Definition: elements.h:3174
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
Definition: elements.cc:4440
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
//////////////////////////////////////////////////////////////////////// ////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
static double two_norm(const Vector< double > &v)
2-norm of a vector
static void scalar_times_vector(const double &lambda, const Vector< double > &v, Vector< double > &lambda_times_v)
Multiply a vector by a scalar.
static double scalar_product(const Vector< double > &v1, const Vector< double > &v2)
Scalar product between two vectors.
double interpolated_u(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
unsigned self_test()
Self-test: Return 0 for OK.
static void allocate_vector_of_vectors(unsigned n_rows, unsigned n_cols, Vector< Vector< double >> &v)
Helper fct: Allocate storage for a vector of vectors of doubles to v(n_rows,n_cols) and initialise ea...
virtual int u_local_eqn(const unsigned &n)
Get the local equation number of the (one and only) unknown stored at local node n (returns -1 if val...
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln at n_plot^2 plot points.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void output(std::ostream &outfile)
Output with default number of plot points.
void exact_position(const Vector< double > &s, Vector< double > &r, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Get exact position vector to meniscus at local coordinate s.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Add the element's contribution to its residual vector.
virtual void get_spine_base(const Vector< double > &x, Vector< double > &spine_base, Vector< Vector< double >> &dspine_base) const
Get spine base vector field: Defaults to standard cartesian representation if no spine base fct point...
virtual double u(const unsigned &n) const
Access function: Nodal function value at local node n Uses suitably interpolated value for hanging no...
double get_kappa() const
Get curvature.
void get_spine(const Vector< double > &x, Vector< double > &spine, Vector< Vector< double >> &dspine) const
Get spine vector field: Defaults to standard cartesian representation if no spine base fct pointers h...
bool use_spines() const
Use spines or not? (Based on availability of function pointers to to spine and spine base vector fiel...
void position(const Vector< double > &s, Vector< double > &r) const
Get position vector to meniscus at local coordinate s.
static void cross_product(const Vector< double > &v1, const Vector< double > &v2, Vector< double > &v_cross)
Cross-product: v_cross= v1 x v2.
static void vector_sum(const Vector< double > &v1, const Vector< double > &v2, Vector< double > &vs)
Vectorial sum of two vectors.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...