~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to Mesh/Field.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-02 18:12:15 UTC
  • mfrom: (1.2.8 upstream)
  • Revision ID: james.westby@ubuntu.com-20090902181215-yla8zvcas2ucvkm9
[Christophe Prud'homme]
* New upstream release
  + fixed surface mesh orientation bug introduced in 2.4.0;
  + mesh and graphics code refactoring;
  + small usability enhancements and bug fixes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
7
7
//   Jonathan Lambrechts
8
8
//
9
9
 
 
10
#include <stdlib.h>
10
11
#include <list>
11
12
#include <math.h>
12
13
#include <fstream>
14
15
#include <string.h>
15
16
#include <sstream>
16
17
#include "GmshConfig.h"
17
 
 
18
 
#if defined(HAVE_MATH_EVAL)
19
 
#include "matheval.h"
20
 
#endif
21
 
#if defined(HAVE_ANN)
22
 
#include "ANN/ANN.h"
23
 
#endif
24
 
 
25
18
#include "Context.h"
26
19
#include "Field.h"
27
20
#include "GeoInterpolation.h"
28
21
#include "GModel.h"
29
22
#include "GmshMessage.h"
30
23
#include "Numeric.h"
 
24
#include "mathEvaluator.h"
31
25
 
32
26
#if !defined(HAVE_NO_POST)
33
27
#include "OctreePost.h"
35
29
#include "MVertex.h"
36
30
#endif
37
31
 
 
32
#if defined(HAVE_ANN)
 
33
#include "ANN/ANN.h"
 
34
#endif
 
35
 
38
36
class FieldOptionDouble : public FieldOption
39
37
{
40
38
 public:
462
460
    options["ZMin"] = new FieldOptionDouble
463
461
      (z_min, "Minimum Z coordinate of the box");
464
462
    options["ZMax"] = new FieldOptionDouble
465
 
      (z_max, "Maximum Z coordinate of the box");
 
463
      (z_max, "Maximum Z coordinate of the box");    
466
464
  }
467
465
  const char *getName()
468
466
  {
469
467
    return "Box";
470
468
  }
471
469
  double operator() (double x, double y, double z, GEntity *ge=0)
472
 
  {
 
470
  {    
473
471
    return (x <= x_max && x >= x_min && y <= y_max && y >= y_min && z <= z_max
474
472
            && z >= z_min) ? v_in : v_out;
475
473
  }
476
474
};
477
475
 
 
476
class CylinderField : public Field
 
477
{
 
478
  double v_in, v_out;
 
479
  double xc,yc,zc;
 
480
  double xa,ya,za;
 
481
  double R;
 
482
  
 
483
 public:
 
484
  std::string getDescription()
 
485
  {
 
486
    return "The value of this field is VIn inside a frustrated cylinder, VOut outside.\n"
 
487
      "The cylinder is given by\n\n"
 
488
      "  ||dX||^2 < R^2 &&\n"
 
489
      "  (X-X0).A < ||A||^2\n"
 
490
      "  dX = (X - X0) - ((X - X0).A)/(||A||^2) . A";
 
491
  }
 
492
  CylinderField()
 
493
  {
 
494
    v_in = v_out = xc = yc = zc = xa = ya = R = 0;
 
495
    za = 1.;
 
496
    
 
497
    options["VIn"] = new FieldOptionDouble
 
498
      (v_in, "Value inside the cylinder");
 
499
    options["VOut"] = new FieldOptionDouble
 
500
      (v_out, "Value outside the cylinder");
 
501
 
 
502
    options["XCenter"] = new FieldOptionDouble
 
503
      (xc, "X coordinate of the cylinder center");
 
504
    options["YCenter"] = new FieldOptionDouble
 
505
      (yc, "Y coordinate of the cylinder center");
 
506
    options["ZCenter"] = new FieldOptionDouble
 
507
      (zc, "Z coordinate of the cylinder center");
 
508
 
 
509
    
 
510
    options["XAxis"] = new FieldOptionDouble
 
511
      (xa, "X component of the cylinder axis");
 
512
    options["YAxis"] = new FieldOptionDouble
 
513
      (ya, "Y component of the cylinder axis");
 
514
    options["ZAxis"] = new FieldOptionDouble
 
515
      (za, "Z component of the cylinder axis");
 
516
 
 
517
    options["Radius"] = new FieldOptionDouble
 
518
      (R,"Radius");    
 
519
    
 
520
  }
 
521
  const char *getName()
 
522
  {
 
523
    return "Cylinder";
 
524
  }
 
525
  double operator() (double x, double y, double z, GEntity *ge=0)
 
526
  {    
 
527
    double dx = x-xc;
 
528
    double dy = y-yc;
 
529
    double dz = z-zc;
 
530
    
 
531
    double adx = (xa * dx + ya * dy + za * dz)/(xa*xa + ya*ya + za*za);
 
532
    
 
533
    dx -= adx * xa;
 
534
    dy -= adx * ya;
 
535
    dz -= adx * za;
 
536
    
 
537
    return ((dx*dx + dy*dy + dz*dz < R*R) && fabs(adx) < 1) ? v_in : v_out;
 
538
    
 
539
  }
 
540
};
 
541
 
478
542
class ThresholdField : public Field
479
543
{
 
544
 protected :
480
545
  int iField;
481
546
  double dmin, dmax, lcmin, lcmax;
482
547
  bool sigmoid, stopAtDistMax;
483
548
 public:
484
 
  const char *getName()
 
549
  virtual const char *getName()
485
550
  {
486
551
    return "Threshold";
487
552
  }
488
 
  std::string getDescription()
 
553
  virtual std::string getDescription()
489
554
  {
490
555
    return "F = LCMin if Field[IField] <= DistMin,\n"
491
556
      "F = LCMax if Field[IField] >= DistMax,\n"
538
603
  }
539
604
};
540
605
 
 
606
class BoundaryLayerField : public ThresholdField {
 
607
public:
 
608
  BoundaryLayerField() 
 
609
  {  }
 
610
  virtual bool isotropic () const {return false;}
 
611
  virtual const char *getName()
 
612
  {
 
613
    return "BoundaryLayer";
 
614
  }
 
615
  virtual std::string getDescription()
 
616
  {
 
617
    return "F = LCMin if Field[IField] <= DistMin,\n"
 
618
      "F = LCMax if Field[IField] >= DistMax,\n"
 
619
      "F = interpolation between LcMin and LcMax if DistMin < Field[IField] < DistMax";
 
620
  }
 
621
  virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
 
622
  {
 
623
    Field *field = GModel::current()->getFields()->get(iField);
 
624
    if(!field) {
 
625
      metr(0,0) = 1/(MAX_LC*MAX_LC);
 
626
      metr(1,1) = 1/(MAX_LC*MAX_LC);
 
627
      metr(2,2) = 1/(MAX_LC*MAX_LC);
 
628
      metr(0,1) = metr(0,2) = metr(1,2) = 0;
 
629
      return;
 
630
    }
 
631
    double dist = (*field) (x, y, z);
 
632
 
 
633
    double r = (dist - dmin) / (dmax - dmin);
 
634
    r = std::max(std::min(r, 1.), 0.);
 
635
    double lc;
 
636
    if(stopAtDistMax && r >= 1.){
 
637
      lc = MAX_LC;
 
638
    }
 
639
    else if(sigmoid){
 
640
      double s = exp(12. * r - 6.) / (1. + exp(12. * r - 6.));
 
641
      lc = lcmin * (1. - s) + lcmax * s;
 
642
    }
 
643
    else{ // linear
 
644
      lc = lcmin * (1 - r) + lcmax * r;
 
645
    }
 
646
    
 
647
    double delta = std::min(CTX::instance()->lc / 1e4, dist);
 
648
    double gx =
 
649
      ((*field) (x + delta / 2, y, z) -
 
650
       (*field) (x - delta / 2, y, z)) / delta;
 
651
    double gy =
 
652
      ((*field) (x, y + delta / 2, z) -
 
653
       (*field) (x, y - delta / 2, z)) / delta;
 
654
    double gz =
 
655
      ((*field) (x, y, z + delta / 2) -
 
656
       (*field) (x, y, z - delta / 2)) / delta;
 
657
 
 
658
    SVector3 g(gx,gy,gz);
 
659
    g.normalize();
 
660
    SVector3 t1,t2;
 
661
 
 
662
    if (fabs(gx) < fabs(gy) && fabs(gx) < fabs(gz))
 
663
      t1 = SVector3(1,0,0);
 
664
    else if (fabs(gy) < fabs(gx) && fabs(gy) < fabs(gz))
 
665
      t1 = SVector3(0,1,0);
 
666
    else
 
667
      t1 = SVector3(0,0,1);
 
668
    
 
669
    t2 = crossprod(g,t1);
 
670
    t2.normalize();
 
671
    t1 = crossprod(t2,g);
 
672
    
 
673
    metr  = SMetric3(1./(lc*lc), 
 
674
                     1/(lcmax*lcmax),
 
675
                     1/(lcmax*lcmax),
 
676
                     g,t1,t2);
 
677
  }
 
678
};
 
679
 
541
680
class GradientField : public Field
542
681
{
543
682
  int iField, kind;
594
733
      return sqrt(gx * gx + gy * gy + gz * gz);
595
734
    default:
596
735
      Msg::Error("Field %i : Unknown kind (%i) of gradient.", this->id,
597
 
                 kind);
 
736
                 kind);
598
737
      return MAX_LC;
599
738
    }
600
739
  }
601
740
};
602
741
 
 
742
 
603
743
class CurvatureField : public Field
604
744
{
605
745
  int iField;
645
785
    grad_norm(*field, x, y, z + delta / 2, grad[4]);
646
786
    grad_norm(*field, x, y, z - delta / 2, grad[5]);
647
787
    return (grad[0][0] - grad[1][0] + grad[2][1] - 
648
 
            grad[3][1] + grad[4][2] - grad[5][2]) / delta;
 
788
            grad[3][1] + grad[4][2] - grad[5][2]) / delta;
649
789
  }
650
790
};
651
791
 
730
870
    Field *field = GModel::current()->getFields()->get(iField);
731
871
    if(!field) return MAX_LC;
732
872
    return ((*field) (x + delta , y, z)+ (*field) (x - delta , y, z)
733
 
            +(*field) (x, y + delta , z)+ (*field) (x, y - delta , z)
734
 
            +(*field) (x, y, z + delta )+ (*field) (x, y, z - delta )
735
 
            -6* (*field) (x , y, z)) / (delta*delta);
 
873
            +(*field) (x, y + delta , z)+ (*field) (x, y - delta , z)
 
874
            +(*field) (x, y, z + delta )+ (*field) (x, y, z - delta )
 
875
            -6* (*field) (x , y, z)) / (delta*delta);
736
876
  }
737
877
};
738
878
 
766
906
    Field *field = GModel::current()->getFields()->get(iField);
767
907
    if(!field) return MAX_LC;
768
908
    return ((*field) (x + delta , y, z) + (*field) (x - delta, y, z)
769
 
            + (*field) (x, y + delta, z) + (*field) (x, y - delta, z)
770
 
            + (*field) (x, y, z + delta) + (*field) (x, y, z - delta)
771
 
            + (*field) (x, y, z)) / 7;
 
909
            + (*field) (x, y + delta, z) + (*field) (x, y - delta, z)
 
910
            + (*field) (x, y, z + delta) + (*field) (x, y, z - delta)
 
911
            + (*field) (x, y, z)) / 7;
772
912
  }
773
913
};
774
914
 
775
 
#if defined(HAVE_MATH_EVAL)
776
915
class MathEvalExpression
777
916
{
778
 
  bool error_status;
779
 
  std::list<Field*> *list;
780
 
  int nvalues;
781
 
  char **names;
782
 
  double *values;
783
 
  void *eval;
784
 
  int *evaluators_id;
785
 
  std::string function;
786
 
  char *c_str_function;
 
917
 private:
 
918
  mathEvaluator *_f;
 
919
  std::set<int> _fields;
787
920
 public:
788
 
  double evaluate(double x, double y, double z)
 
921
  MathEvalExpression() : _f(0) {}
 
922
  ~MathEvalExpression(){ if(_f) delete _f; }
 
923
  bool set_function(const std::string &f)
789
924
  {
790
 
    if(error_status)
791
 
      return MAX_LC;
792
 
    for(int i = 0; i < nvalues; i++){
793
 
      switch (evaluators_id[i]) {
794
 
      case -1:
795
 
        values[i] = x;
796
 
        break;
797
 
      case -2: 
798
 
        values[i] = y;
799
 
        break;
800
 
      case -3: 
801
 
        values[i] = z;
802
 
        break;
803
 
      default:
804
 
        {
805
 
          Field *f = GModel::current()->getFields()->get(evaluators_id[i]);
806
 
          values[i] = f ? (*f) (x, y, z) : MAX_LC;
 
925
    // get id numbers of fields appearing in the function
 
926
    _fields.clear();
 
927
    unsigned int i = 0;
 
928
    while(i < f.size()){
 
929
      unsigned int j = 0;
 
930
      if(f[i] == 'F'){
 
931
        std::string id("");
 
932
        while(i + 1 + j < f.size() && f[i + 1 + j] >= '0' && f[i + 1 + j] <= '9'){
 
933
          id += f[i + 1 + j];
 
934
          j++;
807
935
        }
 
936
        _fields.insert(atoi(id.c_str()));
808
937
      }
809
 
    }
810
 
    return evaluator_evaluate(eval, nvalues, names, values);
811
 
  }
812
 
  MathEvalExpression()
813
 
  {
814
 
    eval = 0;
815
 
    values = 0;
816
 
    c_str_function = 0;
817
 
    evaluators_id = 0;
818
 
  }
819
 
  bool set_function(const std::string & f)
820
 
  {
821
 
    free_members();
822
 
    error_status = false;
823
 
    c_str_function = strdup(f.c_str());
824
 
    eval = evaluator_create(c_str_function);
825
 
    if(!eval) {
826
 
      error_status = true;
 
938
      i += j + 1;
 
939
    }
 
940
    std::vector<std::string> expressions(1), variables(3 + _fields.size());
 
941
    expressions[0] = f;
 
942
    variables[0] = "x"; 
 
943
    variables[1] = "y"; 
 
944
    variables[2] = "z";
 
945
    i = 3;
 
946
    for(std::set<int>::iterator it = _fields.begin(); it != _fields.end(); it++){
 
947
      std::ostringstream sstream;
 
948
      sstream << "F" << *it;
 
949
      variables[i++] = sstream.str();
 
950
    }
 
951
    if(_f) delete _f;
 
952
    _f = new mathEvaluator(expressions, variables);
 
953
    if(expressions.empty()) {
 
954
      delete _f;
 
955
      _f = 0;
827
956
      return false;
828
957
    }
829
 
    evaluator_get_variables(eval, &names, &nvalues);
830
 
    values = new double[nvalues];
831
 
    evaluators_id = new int[nvalues];
832
 
    for(int i = 0; i < nvalues; i++) {
833
 
      int id;
834
 
      if(!strcmp("x", names[i]))
835
 
        evaluators_id[i] = -1;
836
 
      else if(!strcmp("y", names[i]))
837
 
        evaluators_id[i] = -2;
838
 
      else if(!strcmp("z", names[i]))
839
 
        evaluators_id[i] = -3;
840
 
      else if(sscanf(names[i], "F%i", &id) == 1)
841
 
        evaluators_id[i] = id;
842
 
      else {
843
 
        Msg::Error("Unknown matheval argument \"%s\"\n", names[i]);
844
 
        error_status = true;
845
 
        return false;
846
 
      }
847
 
    }
848
958
    return true;
849
959
  }
850
 
  void free_members()
851
 
  {
852
 
    if(c_str_function)
853
 
      free(c_str_function);
854
 
    if(eval)
855
 
      evaluator_destroy(eval);
856
 
    if(values)
857
 
      delete[]values;
858
 
    if(evaluators_id)
859
 
      delete evaluators_id;
860
 
  }
861
 
  ~MathEvalExpression()
862
 
  {
863
 
    free_members();
 
960
  double evaluate(double x, double y, double z)
 
961
  {
 
962
    if(!_f) return MAX_LC;
 
963
    std::vector<double> values(3 + _fields.size()), res(1);
 
964
    values[0] = x;
 
965
    values[1] = y;
 
966
    values[2] = z;
 
967
    int i = 3;
 
968
    for(std::set<int>::iterator it = _fields.begin(); it != _fields.end(); it++){
 
969
      Field *field = GModel::current()->getFields()->get(*it);
 
970
      values[i++] = field ? (*field)(x, y, z) : MAX_LC;
 
971
    }
 
972
    if(_f->eval(values, res)) 
 
973
      return res[0];
 
974
    else
 
975
      return MAX_LC;
864
976
  }
865
977
};
866
978
 
880
992
    if(update_needed) {
881
993
      if(!expr.set_function(f))
882
994
        Msg::Error("Field %i: Invalid matheval expression \"%s\"",
883
 
                   this->id, f.c_str());
 
995
                   this->id, f.c_str());
884
996
      update_needed = false;
885
997
    }
886
998
    return expr.evaluate(x, y, z);
928
1040
      for(int i = 0; i < 3; i++) {
929
1041
        if(!expr[i].set_function(f[i]))
930
1042
          Msg::Error("Field %i : Invalid matheval expression \"%s\"",
931
 
                     this->id, f[i].c_str());
 
1043
                     this->id, f[i].c_str());
932
1044
      }
933
1045
      update_needed = false;
934
1046
    }
935
1047
    Field *field = GModel::current()->getFields()->get(ifield);
936
1048
    if(!field) return MAX_LC;
937
1049
    return (*field)(expr[0].evaluate(x, y, z),
938
 
                    expr[1].evaluate(x, y, z),
939
 
                    expr[2].evaluate(x, y, z));
 
1050
                    expr[1].evaluate(x, y, z),
 
1051
                    expr[2].evaluate(x, y, z));
940
1052
  }
941
1053
  const char *getName()
942
1054
  {
943
1055
    return "Param";
944
1056
  }
945
1057
};
946
 
#endif
947
1058
 
948
1059
#if !defined(HAVE_NO_POST)
949
1060
class PostViewField : public Field
1002
1113
  MinField()
1003
1114
  {
1004
1115
    options["FieldsList"] = new FieldOptionList
1005
 
      (idlist, "Field indices", &update_needed);
 
1116
      (idlist, "Field indices", &update_needed);
1006
1117
  }
1007
1118
  std::string getDescription()
1008
1119
  {
1142
1253
        delete kdtree;
1143
1254
      }
1144
1255
      int totpoints = nodes_id.size() + n_nodes_by_edge * edges_id.size() + 
1145
 
        n_nodes_by_edge * n_nodes_by_edge * faces_id.size();
 
1256
        n_nodes_by_edge * n_nodes_by_edge * faces_id.size();
1146
1257
      if(totpoints)
1147
1258
        zeronodes = annAllocPts(totpoints, 4);
1148
1259
      int k = 0;
1198
1309
        Surface *s = FindSurface(*it);
1199
1310
        if(s) {
1200
1311
          for(int i = 0; i < n_nodes_by_edge; i++) {
1201
 
            for(int j = 0; j < n_nodes_by_edge; j++) {
1202
 
              double u = (double)i / (n_nodes_by_edge - 1);
1203
 
              double v = (double)j / (n_nodes_by_edge - 1);
1204
 
              Vertex V = InterpolateSurface(s, u, v, 0, 0);
1205
 
              zeronodes[k][0] = V.Pos.X;
1206
 
              zeronodes[k][1] = V.Pos.Y;
1207
 
              zeronodes[k++][2] = V.Pos.Z;
1208
 
            }
1209
 
          }
 
1312
            for(int j = 0; j < n_nodes_by_edge; j++) {
 
1313
              double u = (double)i / (n_nodes_by_edge - 1);
 
1314
              double v = (double)j / (n_nodes_by_edge - 1);
 
1315
              Vertex V = InterpolateSurface(s, u, v, 0, 0);
 
1316
              zeronodes[k][0] = V.Pos.X;
 
1317
              zeronodes[k][1] = V.Pos.Y;
 
1318
              zeronodes[k++][2] = V.Pos.Z;
 
1319
            }
 
1320
          }
1210
1321
        }
1211
1322
        else {
1212
1323
          GFace *f = GModel::current()->getFaceByTag(*it);
1213
1324
          if(f) {
1214
1325
            for(int i = 0; i < n_nodes_by_edge; i++) {
1215
 
              for(int j = 0; j < n_nodes_by_edge; j++) {
1216
 
                double u = (double)i / (n_nodes_by_edge - 1);
1217
 
                double v = (double)j / (n_nodes_by_edge - 1);
1218
 
                Range<double> b1 = ge->parBounds(0);
1219
 
                Range<double> b2 = ge->parBounds(1);
1220
 
                double t1 = b1.low() + u * (b1.high() - b1.low());
1221
 
                double t2 = b2.low() + v * (b2.high() - b2.low());
1222
 
                GPoint gp = f->point(t1, t2);
1223
 
                zeronodes[k][0] = gp.x();
1224
 
                zeronodes[k][1] = gp.y();
1225
 
                zeronodes[k++][2] = gp.z();
1226
 
              }
1227
 
            }
 
1326
              for(int j = 0; j < n_nodes_by_edge; j++) {
 
1327
                double u = (double)i / (n_nodes_by_edge - 1);
 
1328
                double v = (double)j / (n_nodes_by_edge - 1);
 
1329
                Range<double> b1 = ge->parBounds(0);
 
1330
                Range<double> b2 = ge->parBounds(1);
 
1331
                double t1 = b1.low() + u * (b1.high() - b1.low());
 
1332
                double t2 = b2.low() + v * (b2.high() - b2.low());
 
1333
                GPoint gp = f->point(t1, t2);
 
1334
                zeronodes[k][0] = gp.x();
 
1335
                zeronodes[k][1] = gp.y();
 
1336
                zeronodes[k++][2] = gp.z();
 
1337
              }
 
1338
            }
1228
1339
          }
1229
1340
        }
1230
1341
      }
1236
1347
    return sqrt(dist[0]);
1237
1348
  }
1238
1349
};
 
1350
 
1239
1351
#endif
1240
1352
 
1241
1353
template<class F> class FieldFactoryT : public FieldFactory {
1252
1364
{
1253
1365
  map_type_name["Structured"] = new FieldFactoryT<StructuredField>();
1254
1366
  map_type_name["Threshold"] = new FieldFactoryT<ThresholdField>();
 
1367
  map_type_name["BoundaryLayer"] = new FieldFactoryT<BoundaryLayerField>();
1255
1368
  map_type_name["Box"] = new FieldFactoryT<BoxField>();
 
1369
  map_type_name["Cylinder"] = new FieldFactoryT<CylinderField>();
1256
1370
  map_type_name["LonLat"] = new FieldFactoryT<LonLatField>();
1257
1371
#if !defined(HAVE_NO_POST)
1258
1372
  map_type_name["PostView"] = new FieldFactoryT<PostViewField>();
1265
1379
  map_type_name["Laplacian"] = new FieldFactoryT<LaplacianField>();
1266
1380
  map_type_name["Mean"] = new FieldFactoryT<MeanField>();
1267
1381
  map_type_name["Curvature"] = new FieldFactoryT<CurvatureField>();
1268
 
#if defined(HAVE_MATH_EVAL)
1269
1382
  map_type_name["Param"] = new FieldFactoryT<ParametricField>();
1270
1383
  map_type_name["MathEval"] = new FieldFactoryT<MathEvalField>();
1271
 
#endif
1272
1384
#if defined(HAVE_ANN)
1273
1385
  map_type_name["Attractor"] = new FieldFactoryT<AttractorField>();
1274
1386
#endif
1309
1421
    for(int ele = 0; ele < data->getNumElements(0, ent); ele++){
1310
1422
      if(data->skipElement(0, ent, ele)) continue;
1311
1423
      for(int nod = 0; nod < data->getNumNodes(0, ent, ele); nod++){
1312
 
        double x, y, z;
1313
 
        data->getNode(0, ent, ele, nod, x, y, z);
 
1424
        double x, y, z;
 
1425
        data->getNode(0, ent, ele, nod, x, y, z);
1314
1426
        double val = (*this)(x, y, z);
1315
 
        for(int comp = 0; comp < data->getNumComponents(0, ent, ele); comp++)
1316
 
          data->setValue(0, ent, ele, nod, comp, val);
 
1427
        for(int comp = 0; comp < data->getNumComponents(0, ent, ele); comp++)
 
1428
          data->setValue(0, ent, ele, nod, comp, val);
1317
1429
      }
1318
1430
    }
1319
1431
  }