← Back to branch summary

~chaffra/+junk/trilinos

« back to all changes in this revision

Viewing changes to packages/Sundance/tests-std-framework/Problem/Poisson2D.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Johannes Ring
  • Date: 2009-12-13 12:53:22 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20091213125322-in0nrdjc55deqsw9
Tags: 10.0.3.dfsg-1
[Christophe Prud'homme]
* New upstream release

[Johannes Ring]
* debian/patches/libname.patch: Add prefix 'libtrilinos_' to all
  libraries. 
* debian/patches/soname.patch: Add soversion to libraries.
* debian/watch: Update download URL.
* debian/control:
  - Remove python-numeric from Build-Depends (virtual package).
  - Remove automake and autotools from Build-Depends and add cmake to
    reflect switch to CMake.
  - Add python-support to Build-Depends.
* debian/rules: 
  - Cleanup and updates for switch to CMake.

Show diffs side-by-side

added added

removed removed

Lines of Context:
packages/Sundance/tests-std-framework/Problem/Poisson2D.cpp
1
 
/* @HEADER@ */
2
 
// ************************************************************************
3
 
// 
4
 
//                              Sundance
5
 
//                 Copyright (2005) Sandia Corporation
6
 
// 
7
 
// Copyright (year first published) Sandia Corporation.  Under the terms 
8
 
// of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
9
 
// retains certain rights in this software.
10
 
// 
11
 
// This library is free software; you can redistribute it and/or modify
12
 
// it under the terms of the GNU Lesser General Public License as
13
 
// published by the Free Software Foundation; either version 2.1 of the
14
 
// License, or (at your option) any later version.
15
 
//  
16
 
// This library is distributed in the hope that it will be useful, but
17
 
// WITHOUT ANY WARRANTY; without even the implied warranty of
18
 
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19
 
// Lesser General Public License for more details.
20
 
//                                                                                 
21
 
// You should have received a copy of the GNU Lesser General Public
22
 
// License along with this library; if not, write to the Free Software
23
 
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
24
 
// USA                                                                                
25
 
// Questions? Contact Kevin Long (krlong@sandia.gov), 
26
 
// Sandia National Laboratories, Livermore, California, USA
27
 
// 
28
 
// ************************************************************************
29
 
/* @HEADER@ */
 
1
/* 
 
2
 * <Ignore> 
 
3
 * 
 
4
 *</Ignore>
 
5
 */
30
6
 
 
7
/* <Ignore> */
31
8
#include "Sundance.hpp"
32
 
#include "SundanceEvaluator.hpp"
33
 
#include "SundanceExodusMeshReader.hpp"
34
 
 
35
9
using SundanceCore::List;
36
 
/** 
 
10
 
 
11
#if defined(HAVE_SUNDANCE_EXODUS)
 
12
 
 
13
/* </Ignore> */
 
14
 
 
15
/*
 
16
 * \documentclass[10pt]{article}
 
17
 * \begin{document}
 
18
 * 
 
19
 */
 
20
 
 
21
 
 
22
/* 
37
23
 * Solves the Poisson equation in 2D
38
24
 */
39
25
 
42
28
CELL_PREDICATE(RightPointTest, {return fabs(x[0]-1.0) < 1.0e-10;})
43
29
CELL_PREDICATE(TopPointTest, {return fabs(x[1]-1.0) < 1.0e-10;}) 
44
30
 
45
 
#ifdef HAVE_EXODUS
46
31
 
47
32
 
48
33
int main(int argc, char** argv)
49
34
{
50
 
  
51
35
  try
52
 
                {
53
 
      int nx = 2;
54
 
      int ny = 2;
55
 
      string meshFile="builtin";
56
 
      string solverFile = "aztec-ml.xml";
57
 
      Sundance::setOption("meshFile", meshFile, "mesh file");
58
 
      Sundance::setOption("nx", nx, "number of elements in x");
59
 
      Sundance::setOption("ny", ny, "number of elements in y");
60
 
      Sundance::setOption("solver", solverFile, "name of XML file for solver");
61
 
 
62
 
      Sundance::init(&argc, &argv);
63
 
      int np = MPIComm::world().getNProc();
64
 
 
65
 
      /* We will do our linear algebra using Epetra */
66
 
      VectorType<double> vecType = new EpetraVectorType();
67
 
 
68
 
      /* Create a mesh. It will be of type BasisSimplicialMesh, and will
69
 
       * be built using a PartitionedRectangleMesher. */
70
 
      MeshType meshType = new BasicSimplicialMeshType();
71
 
      
72
 
      MeshSource mesher;
73
 
      if (meshFile != "builtin")
74
 
      {
75
 
        mesher = new ExodusMeshReader("../../../examples-tutorial/meshes/"+meshFile, meshType);
76
 
      }
77
 
      else
78
 
      {
79
 
        mesher = new PartitionedRectangleMesher(0.0, 1.0, nx, np, 
80
 
          0.0,  1.0, ny, 1, meshType);
81
 
      }
82
 
      Mesh mesh = mesher.getMesh();
83
 
 
84
 
 
85
 
      bool meshOK = mesh.checkConsistency(meshFile+"-check");
86
 
      if (meshOK) 
87
 
      {
88
 
        cout << "mesh is OK" << endl;
89
 
      }
90
 
      else
91
 
      {
92
 
        cout << "mesh is INCONSISTENT" << endl;
93
 
      }
94
 
      mesh.dump(meshFile+"-dump");
95
 
 
96
 
      /* Create a cell filter that will identify the maximal cells
97
 
       * in the interior of the domain */
98
 
      CellFilter interior = new MaximalCellFilter();
99
 
      CellFilter edges = new DimensionalCellFilter(1);
100
 
 
101
 
      CellFilter left;
102
 
      CellFilter right;
103
 
      CellFilter top;
104
 
      CellFilter bottom;
105
 
 
106
 
      if (meshFile != "builtin")
107
 
      {
108
 
        left = edges.labeledSubset(1);
109
 
        right = edges.labeledSubset(2);
110
 
        top = edges.labeledSubset(3);
111
 
        bottom = edges.labeledSubset(4);
112
 
      }
113
 
      else
114
 
      {
115
 
        left = edges.subset(new LeftPointTest());
116
 
        right = edges.subset(new RightPointTest());
117
 
        top = edges.subset(new TopPointTest());
118
 
        bottom = edges.subset(new BottomPointTest());
119
 
      }
120
 
      
121
 
      /* Create unknown and test functions, discretized using second-order
122
 
       * Lagrange interpolants */
123
 
      BasisFamily basis = new Lagrange(1);
124
 
      Expr u = new UnknownFunction(basis, "u");
125
 
      Expr v = new TestFunction(basis, "v");
126
 
 
127
 
      /* Create differential operator and coordinate functions */
128
 
      Expr dx = new Derivative(0);
129
 
      Expr dy = new Derivative(1);
130
 
      Expr grad = List(dx, dy);
131
 
      Expr x = new CoordExpr(0);
132
 
      Expr y = new CoordExpr(1);
133
 
 
134
 
      /* We need a quadrature rule for doing the integrations */
135
 
      QuadratureFamily quad2 = new GaussianQuadrature(2);
136
 
      QuadratureFamily quad4 = new GaussianQuadrature(4);
137
 
 
138
 
      /* Define the weak form */
139
 
      //Expr eqn = Integral(interior, (grad*v)*(grad*u) + v, quad);
140
 
      Expr one = new SundanceCore::Parameter(1.0);
141
 
      Expr exactSoln = x+2.0*y;//0.5*x*x + (1.0/3.0)*y;
142
 
      Expr eqn = Integral(interior, (grad*u)*(grad*v)  /* + one*v */, quad2);
143
 
      /* Define the Dirichlet BC */
144
 
      Expr h = new CellDiameterExpr();
145
 
      Expr bc = EssentialBC(bottom+top+left+right, v*(u-exactSoln)/h, quad4);
146
 
 
147
 
      Assembler::defaultVerbParams()->set<int>("global", 2);
148
 
      Assembler::defaultVerbParams()->set<int>("evaluation", 4);
149
 
      Assembler::defaultVerbParams()->set<int>("assembly loop", 2);
150
 
      
151
 
      /* We can now set up the linear problem! */
152
 
      LinearProblem prob(mesh, eqn, bc, v, u, vecType);
153
 
 
154
 
      ParameterXMLFileReader reader(searchForFile("SolverParameters/" + solverFile));
155
 
      ParameterList solverParams = reader.getParameters();
156
 
      LinearSolver<double> solver 
157
 
        = LinearSolverBuilder::createSolver(solverParams);
158
 
 
159
 
      //      cout << "map = " << endl;
160
 
      //prob.rowMap()->print(cout);
161
 
      //cout << endl;
162
 
      Expr soln = prob.solve(solver);
163
 
 
164
 
      DiscreteSpace discSpace2(mesh, new Lagrange(2), vecType);
165
 
      DiscreteSpace discSpace0(mesh, new Lagrange(0), vecType);
166
 
      L2Projector proj1(discSpace2, soln-exactSoln);
167
 
      double pid = MPIComm::world().getRank();
168
 
      L2Projector proj2(discSpace0, pid);
169
 
      Expr errorDisc = proj1.project();
170
 
      Expr pidDisc = proj2.project();
171
 
 
172
 
 
173
 
 
174
 
 
175
 
      /* Write the field in VTK format */
176
 
      FieldWriter w = new VTKWriter("Poisson2d");
177
 
      w.addMesh(mesh);
178
 
      w.addField("soln", new ExprFieldWrapper(soln[0]));
179
 
      w.addField("error", new ExprFieldWrapper(errorDisc));
180
 
      w.addField("rank", new ExprFieldWrapper(pidDisc));
181
 
      w.write();
182
 
 
183
 
      Expr err = exactSoln - soln;
184
 
      Expr errExpr = Integral(interior, 
185
 
                              err*err,
186
 
                              quad4);
187
 
 
188
 
      Expr derivErr = dx*(exactSoln-soln);
189
 
      Expr derivErrExpr = Integral(interior, 
190
 
                                   derivErr*derivErr, 
191
 
                                   quad2);
192
 
 
193
 
      
194
 
 
195
 
      Expr fluxErrExpr = Integral(top, 
196
 
                                  pow(dy*(soln-exactSoln), 2),
197
 
                                  new GaussianQuadrature(2));
198
 
 
199
 
      
200
 
 
201
 
      Expr exactFluxExpr = Integral(top, 
202
 
                                    dy*exactSoln,
203
 
                                    new GaussianQuadrature(2));
204
 
 
205
 
      Expr numFluxExpr = Integral(top, 
206
 
                                  dy*soln,
207
 
                                  new GaussianQuadrature(2));
208
 
      //        + Integral(top+bottom, 
209
 
      //                   pow(dy*(soln-exactSoln), 2),
210
 
      //                   new GaussianQuadrature(2));
211
 
 
212
 
 
213
 
      FunctionalEvaluator errInt(mesh, errExpr);
214
 
      FunctionalEvaluator derivErrInt(mesh, derivErrExpr);
215
 
 
216
 
      double errorSq = errInt.evaluate();
217
 
      cout << "error norm = " << sqrt(errorSq) << endl << endl;
218
 
 
219
 
      double derivErrorSq = derivErrInt.evaluate();
220
 
      cout << "deriv error norm = " << sqrt(derivErrorSq) << endl << endl;
221
 
 
222
 
      Assembler::defaultVerbParams()->set<int>("global", 2);
223
 
      Assembler::defaultVerbParams()->set<int>("assembly loop", 2);
224
 
      double fluxErrorSq = evaluateIntegral(mesh, fluxErrExpr);
225
 
      cout << "flux error norm = " << sqrt(fluxErrorSq) << endl << endl;
226
 
 
227
 
      cout << "exact flux = " << evaluateIntegral(mesh, exactFluxExpr) << endl;
228
 
      cout << "numerical flux = " << evaluateIntegral(mesh, numFluxExpr) << endl;
229
 
 
230
 
      Sundance::passFailTest(sqrt(errorSq + derivErrorSq + fluxErrorSq), 1.0e-9);
231
 
 
232
 
    }
 
36
  {
 
37
    int nx = 2;
 
38
    int ny = 2;
 
39
    string meshFile="builtin";
 
40
    string solverFile = "aztec-ml.xml";
 
41
    Sundance::setOption("meshFile", meshFile, "mesh file");
 
42
    Sundance::setOption("nx", nx, "number of elements in x");
 
43
    Sundance::setOption("ny", ny, "number of elements in y");
 
44
    Sundance::setOption("solver", solverFile, "name of XML file for solver");
 
45
 
 
46
    Sundance::init(&argc, &argv);
 
47
    int np = MPIComm::world().getNProc();
 
48
 
 
49
    nx = nx*np;
 
50
    ny = ny*np;
 
51
 
 
52
    /* 
 
53
     * <Header level="subsubsection" name="vector_type">
 
54
     * Creation of vector type
 
55
     * </Header>
 
56
     * 
 
57
     * Next we create a \verb+VectorType+ object. A \verb+VectorType+
 
58
     * is an abstract factory which produces \verb+VectorSpace+ objects.
 
59
     * A \verb+VectorSpace+ is itself a factory which can produce 
 
60
     * \verb+Vector+ objects which are, logically enough, vectors.
 
61
     * 
 
62
     * Why this heirarchy of factorys? Vectors may need to be created
 
63
     * inside code far below the user level, and the \verb+VectorSpace+ 
 
64
     * encapsulates the information needed to build a vector of a 
 
65
     * given type, size, and distribution over processors. The 
 
66
     * \verb+VectorType+ is needed because we might need to select a 
 
67
     * particular {\it software representation} for a vector, {\it e.g.}, 
 
68
     * Trilinos, Petsc, or some other library. 
 
69
     *
 
70
     * By using an \verb+EpetraVectorType+, we will be creating
 
71
     * \verb+EpetraVectorSpace+ vector spaces, which in turn
 
72
     * create \verb+EpetraVector+ vectors.
 
73
     */
 
74
    VectorType<double> vecType = new EpetraVectorType();
 
75
 
 
76
    /* 
 
77
     * <Header level="subsubsection" name="mesh">
 
78
     * Creation of mesh
 
79
     * </Header>
 
80
     * 
 
81
     * The creation of a mesh object involves several intermediate
 
82
     * objects: a \verb+MeshType+ and a \verb+MeshSource+. The 
 
83
     * \verb+MeshType+ specifies what mesh data structure will be used,
 
84
     * and the \verb+MeshSource+ builds that data structure. How it is
 
85
     * built will depend on the \verb+MeshSource+ subtype used. An
 
86
     * \verb+ExodusMeshReader+, for instance, reads from an Exodus II file,
 
87
     * while a \verb+PartitionedRectangleMesher+ builds a uniform 
 
88
     * triangulation of a rectangle.
 
89
     * 
 
90
     * Once the \verb+MeshType+ and \verb+MeshSource+ objects have been
 
91
     * defined, the mesh is created (read, built, whatever) by a call
 
92
     * to the \verb+getMesh()+ method of \verb+MeshSource+.
 
93
     */
 
94
    MeshType meshType = new BasicSimplicialMeshType();
 
95
      
 
96
    MeshSource mesher;
 
97
    if (meshFile != "builtin")
 
98
    {
 
99
      mesher = new ExodusMeshReader("../../../examples-tutorial/meshes/"+meshFile, meshType);
 
100
    }
 
101
    else
 
102
    {
 
103
      int npx = -1;
 
104
      int npy = -1;
 
105
      PartitionedRectangleMesher::balanceXY(np, &npx, &npy);
 
106
      TEST_FOR_EXCEPT(npx < 1);
 
107
      TEST_FOR_EXCEPT(npy < 1);
 
108
      TEST_FOR_EXCEPT(npx * npy != np);
 
109
      mesher = new PartitionedRectangleMesher(0.0, 1.0, nx, npx, 
 
110
        0.0,  1.0, ny, npy, meshType);
 
111
    }
 
112
    Mesh mesh = mesher.getMesh();
 
113
 
 
114
 
 
115
    bool meshOK = mesh.checkConsistency(meshFile+"-check");
 
116
    if (meshOK) 
 
117
    {
 
118
      cout << "mesh is OK" << endl;
 
119
    }
 
120
    else
 
121
    {
 
122
      cout << "mesh is INCONSISTENT" << endl;
 
123
    }
 
124
    mesh.dump(meshFile+"-dump");
 
125
 
 
126
    WatchFlag watchMe("watch eqn");
 
127
    watchMe.deactivate();
 
128
 
 
129
    WatchFlag watchBC("watch BCs");
 
130
    watchBC.setParam("integration setup", 0);
 
131
    watchBC.setParam("integration", 0);
 
132
    watchBC.setParam("fill", 0);
 
133
    watchBC.setParam("evaluation", 0);
 
134
    watchBC.deactivate();
 
135
    
 
136
 
 
137
 
 
138
    /* 
 
139
     * <Header level="subsubsection" name="cell_filter">
 
140
     * Specification of geometric regions
 
141
     * </Header>
 
142
     * 
 
143
     * We'll need to specify subsets of the mesh on which equations
 
144
     * or boundary conditions are defined. In many FEA codes this is
 
145
     * done by explicit definition of element blocks, node sets, and
 
146
     * side sets. Rather than working with sets explicitly at the user
 
147
     * level, we instead work with filtering rules that produce 
 
148
     * sets of cells. These rules are represented by \verb+CellFilter+
 
149
     * objects. 
 
150
     *
 
151
     * \verb+MaximalCellFilter+
 
152
     * selects all cells having dimension equal to the spatial dimension
 
153
     * of the mesh. 
 
154
     */
 
155
    CellFilter interior = new MaximalCellFilter();
 
156
 
 
157
    
 
158
    /* 
 
159
     * \verb+DimensionalCellFilter+
 
160
     * selects all cells of a specified dimension.
 
161
     */
 
162
    CellFilter edges = new DimensionalCellFilter(1);
 
163
 
 
164
    CellFilter left;
 
165
    CellFilter right;
 
166
    CellFilter top;
 
167
    CellFilter bottom;
 
168
 
 
169
    if (meshFile != "builtin")
 
170
    {
 
171
      left = edges.labeledSubset(1);
 
172
      right = edges.labeledSubset(2);
 
173
      top = edges.labeledSubset(3);
 
174
      bottom = edges.labeledSubset(4);
 
175
    }
 
176
    else
 
177
    {
 
178
      left = edges.subset(new LeftPointTest());
 
179
      right = edges.subset(new RightPointTest());
 
180
      top = edges.subset(new TopPointTest());
 
181
      bottom = edges.subset(new BottomPointTest());
 
182
    }
 
183
      
 
184
    /* Create unknown and test functions, discretized using second-order
 
185
     * Lagrange interpolants */
 
186
    BasisFamily basis = new Lagrange(1);
 
187
    Expr u = new UnknownFunction(basis, "u");
 
188
    Expr v = new TestFunction(basis, "v");
 
189
 
 
190
    /* Create differential operator and coordinate functions */
 
191
    Expr dx = new Derivative(0);
 
192
    Expr dy = new Derivative(1);
 
193
    Expr grad = List(dx, dy);
 
194
    Expr x = new CoordExpr(0);
 
195
    Expr y = new CoordExpr(1);
 
196
 
 
197
    /* We need a quadrature rule for doing the integrations */
 
198
    QuadratureFamily quad2 = new GaussianQuadrature(2);
 
199
    QuadratureFamily quad4 = new GaussianQuadrature(4);
 
200
 
 
201
    /* Define the weak form */
 
202
    //Expr eqn = Integral(interior, (grad*v)*(grad*u) + v, quad);
 
203
    Expr one = new SundanceCore::Parameter(1.0);
 
204
    Expr exactSoln = 2.0*x+y;
 
205
    Expr eqn = Integral(interior, (grad*u)*(grad*v), quad2, watchMe);
 
206
    /* Define the Dirichlet BC */
 
207
    Expr h = new CellDiameterExpr();
 
208
    Expr bc = EssentialBC(bottom+top+left+right, v*(u-exactSoln), quad4, watchBC);
 
209
 
 
210
    /* We can now set up the linear problem! */
 
211
    LinearProblem prob(mesh, eqn, bc, v, u, vecType);
 
212
 
 
213
 
 
214
    ParameterXMLFileReader reader(solverFile);
 
215
    ParameterList solverParams = reader.getParameters();
 
216
    LinearSolver<double> solver 
 
217
      = LinearSolverBuilder::createSolver(solverParams);
 
218
 
 
219
    Expr soln = prob.solve(solver);
 
220
 
 
221
    DiscreteSpace discSpace2(mesh, new Lagrange(2), vecType);
 
222
    DiscreteSpace discSpace0(mesh, new Lagrange(0), vecType);
 
223
    L2Projector proj1(discSpace2, soln-exactSoln);
 
224
    double pid = MPIComm::world().getRank();
 
225
    L2Projector proj2(discSpace0, pid);
 
226
    Expr errorDisc = proj1.project();
 
227
    Expr pidDisc = proj2.project();
 
228
 
 
229
 
 
230
 
 
231
 
 
232
    /* Write the field in VTK format */
 
233
    FieldWriter w = new VTKWriter("Poisson2d");
 
234
    w.addMesh(mesh);
 
235
    w.addField("soln", new ExprFieldWrapper(soln[0]));
 
236
    w.addField("error", new ExprFieldWrapper(errorDisc));
 
237
    w.addField("rank", new ExprFieldWrapper(pidDisc));
 
238
    w.write();
 
239
 
 
240
    FieldWriter w2 = new VerboseFieldWriter("mesh");
 
241
    w2.addMesh(mesh);
 
242
    w2.write();
 
243
 
 
244
    Expr err = exactSoln - soln;
 
245
    Expr errExpr = Integral(interior, 
 
246
      err*err,
 
247
      quad4);
 
248
 
 
249
    Expr derivErr = dx*(exactSoln-soln);
 
250
    Expr derivErrExpr = Integral(interior, 
 
251
      derivErr*derivErr, 
 
252
      quad2);
 
253
 
 
254
      
 
255
 
 
256
    Expr fluxErrExpr = Integral(top, 
 
257
      pow(dy*(soln-exactSoln), 2),
 
258
      new GaussianQuadrature(2));
 
259
 
 
260
      
 
261
    watchBC.deactivate();
 
262
    Expr exactFluxExpr = Integral(top, 
 
263
      dy*exactSoln,
 
264
      new GaussianQuadrature(2), watchBC);
 
265
 
 
266
    Expr numFluxExpr = Integral(top, 
 
267
      dy*soln,
 
268
      new GaussianQuadrature(2));
 
269
 
 
270
 
 
271
    FunctionalEvaluator errInt(mesh, errExpr);
 
272
    FunctionalEvaluator derivErrInt(mesh, derivErrExpr);
 
273
 
 
274
    double errorSq = errInt.evaluate();
 
275
    cout << "error norm = " << sqrt(errorSq) << endl << endl;
 
276
 
 
277
    double derivErrorSq = derivErrInt.evaluate();
 
278
    cout << "deriv error norm = " << sqrt(derivErrorSq) << endl << endl;
 
279
 
 
280
 
 
281
    double fluxErrorSq = evaluateIntegral(mesh, fluxErrExpr);
 
282
    cout << "flux error norm = " << sqrt(fluxErrorSq) << endl << endl;
 
283
 
 
284
 
 
285
    cout << "exact flux = " << evaluateIntegral(mesh, exactFluxExpr) << endl;
 
286
    cout << "numerical flux = " << evaluateIntegral(mesh, numFluxExpr) << endl;
 
287
 
 
288
    Sundance::passFailTest(sqrt(errorSq + derivErrorSq + fluxErrorSq), 1.0e-9);
 
289
 
 
290
  }
233
291
        catch(exception& e)
234
 
                {
235
 
      Sundance::handleException(e);
236
 
                }
 
292
  {
 
293
    Sundance::handleException(e);
 
294
  }
237
295
  Sundance::finalize(); return Sundance::testStatus(); 
238
296
 
239
297
  return Sundance::testStatus();
240
298
}
241
299
 
 
300
 
 
301
/* <Ignore> */
242
302
#else
243
303
 
244
304
 
245
305
int main(int argc, char** argv)
246
306
{
 
307
  Sundance::init(&argc, &argv);
247
308
  std::cout << "dummy Poisson2D PASSED. Enable exodus to run the actual test" << std::endl;
 
309
  Sundance::finalize();
248
310
  return 0;
249
311
}
250
312
 
251
313
 
252
314
#endif
 
315
 
 
316
/* </Ignore> */
 
317
 
 
318
/* \end{document} */