rename files curvedelems_new -> curvedelems

bugfix with lock(mesh.mutex) in vsmesh

libsrc/csg/genmesh.cpp

@@ -661,7 +661,6 @@ namespace netgen
 		    int perfstepsstart, int perfstepsend,
 		    const char * optstr)
   {
-
     if (mesh && mesh->GetNSE() &&
 	!geom.GetNSolids())
       {
@@ -673,8 +672,14 @@ namespace netgen
 
     if (perfstepsstart <= MESHCONST_ANALYSE)
       {
+        /*
 	delete mesh;
 	mesh = new Mesh();
+        */
+        if (mesh)
+          mesh -> DeleteMesh();
+        else
+          mesh = new Mesh();
 
 	mesh->SetGlobalH (mparam.maxh);
 	mesh->SetMinimalH (mparam.minh);

libsrc/meshing/Makefile.am

@@ -4,20 +4,20 @@ specials.hpp bisect.hpp geomsearch.hpp hpref_tet.hpp meshing3.hpp     \
 topology.hpp boundarylayer.hpp global.hpp hpref_trig.hpp meshing.hpp  \
 validate.hpp classifyhpel.hpp hpref_hex.hpp improve2.hpp meshtool.hpp \
 clusters.hpp hprefinement.hpp improve3.hpp meshtype.hpp		      \
-curvedelems.hpp hpref_prism.hpp localh.hpp msghandler.hpp	      \
+ hpref_prism.hpp localh.hpp msghandler.hpp curvedelems.hpp	      \
-curvedelems_new.hpp hpref_pyramid.hpp meshclass.hpp ruler2.hpp
+ hpref_pyramid.hpp meshclass.hpp ruler2.hpp
 
 
 AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include
 METASOURCES = AUTO
 noinst_LTLIBRARIES = libmesh.la
-libmesh_la_SOURCES = adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp \
+libmesh_la_SOURCES = adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp    \
-	clusters.cpp curvedelems_new.cpp delaunay.cpp geomsearch.cpp global.cpp \
+	clusters.cpp curvedelems.cpp delaunay.cpp geomsearch.cpp global.cpp    \
-	hprefinement.cpp improve2.cpp improve2gen.cpp improve3.cpp localh.cpp meshclass.cpp \
+	 hprefinement.cpp improve2.cpp improve2gen.cpp improve3.cpp localh.cpp \
-	meshfunc.cpp meshfunc2d.cpp meshing2.cpp meshing3.cpp meshtool.cpp meshtype.cpp \
+	meshclass.cpp meshfunc.cpp meshfunc2d.cpp meshing2.cpp meshing3.cpp    \
-	msghandler.cpp netrule2.cpp netrule3.cpp parser2.cpp parser3.cpp prism2rls.cpp \
+	meshtool.cpp meshtype.cpp msghandler.cpp netrule2.cpp		       \
-	pyramid2rls.cpp pyramidrls.cpp quadrls.cpp refine.cpp ruler2.cpp \
+	netrule3.cpp parser2.cpp parser3.cpp prism2rls.cpp		       \
-	ruler3.cpp secondorder.cpp smoothing2.5.cpp smoothing2.cpp smoothing3.cpp \
+	pyramid2rls.cpp pyramidrls.cpp quadrls.cpp refine.cpp		       \
-	specials.cpp tetrarls.cpp topology.cpp triarls.cpp validate.cpp zrefine.cpp
+	ruler2.cpp ruler3.cpp secondorder.cpp smoothing2.5.cpp		       \
-
+	smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp		       \
-
+	topology.cpp triarls.cpp validate.cpp zrefine.cpp

libsrc/meshing/Makefile.in

@@ -49,7 +49,7 @@ CONFIG_CLEAN_FILES =
 LTLIBRARIES = $(noinst_LTLIBRARIES)
 libmesh_la_LIBADD =
 am_libmesh_la_OBJECTS = adfront2.lo adfront3.lo bisect.lo \
-	boundarylayer.lo clusters.lo curvedelems_new.lo delaunay.lo \
+	boundarylayer.lo clusters.lo curvedelems.lo delaunay.lo \
 	geomsearch.lo global.lo hprefinement.lo improve2.lo \
 	improve2gen.lo improve3.lo localh.lo meshclass.lo meshfunc.lo \
 	meshfunc2d.lo meshing2.lo meshing3.lo meshtool.lo meshtype.lo \
@@ -237,20 +237,22 @@ specials.hpp bisect.hpp geomsearch.hpp hpref_tet.hpp meshing3.hpp     \
 topology.hpp boundarylayer.hpp global.hpp hpref_trig.hpp meshing.hpp  \
 validate.hpp classifyhpel.hpp hpref_hex.hpp improve2.hpp meshtool.hpp \
 clusters.hpp hprefinement.hpp improve3.hpp meshtype.hpp		      \
-curvedelems.hpp hpref_prism.hpp localh.hpp msghandler.hpp	      \
+ hpref_prism.hpp localh.hpp msghandler.hpp curvedelems.hpp	      \
-curvedelems_new.hpp hpref_pyramid.hpp meshclass.hpp ruler2.hpp
+ hpref_pyramid.hpp meshclass.hpp ruler2.hpp
 
 AM_CPPFLAGS = -I$(top_srcdir)/libsrc/include
 METASOURCES = AUTO
 noinst_LTLIBRARIES = libmesh.la
-libmesh_la_SOURCES = adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp \
+libmesh_la_SOURCES = adfront2.cpp adfront3.cpp bisect.cpp boundarylayer.cpp    \
-	clusters.cpp curvedelems_new.cpp delaunay.cpp geomsearch.cpp global.cpp \
+	clusters.cpp curvedelems.cpp delaunay.cpp geomsearch.cpp global.cpp    \
-	hprefinement.cpp improve2.cpp improve2gen.cpp improve3.cpp localh.cpp meshclass.cpp \
+	 hprefinement.cpp improve2.cpp improve2gen.cpp improve3.cpp localh.cpp \
-	meshfunc.cpp meshfunc2d.cpp meshing2.cpp meshing3.cpp meshtool.cpp meshtype.cpp \
+	meshclass.cpp meshfunc.cpp meshfunc2d.cpp meshing2.cpp meshing3.cpp    \
-	msghandler.cpp netrule2.cpp netrule3.cpp parser2.cpp parser3.cpp prism2rls.cpp \
+	meshtool.cpp meshtype.cpp msghandler.cpp netrule2.cpp		       \
-	pyramid2rls.cpp pyramidrls.cpp quadrls.cpp refine.cpp ruler2.cpp \
+	netrule3.cpp parser2.cpp parser3.cpp prism2rls.cpp		       \
-	ruler3.cpp secondorder.cpp smoothing2.5.cpp smoothing2.cpp smoothing3.cpp \
+	pyramid2rls.cpp pyramidrls.cpp quadrls.cpp refine.cpp		       \
-	specials.cpp tetrarls.cpp topology.cpp triarls.cpp validate.cpp zrefine.cpp
+	ruler2.cpp ruler3.cpp secondorder.cpp smoothing2.5.cpp		       \
+	smoothing2.cpp smoothing3.cpp specials.cpp tetrarls.cpp		       \
+	topology.cpp triarls.cpp validate.cpp zrefine.cpp
 
 all: all-am
 
@@ -308,7 +310,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/bisect.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/boundarylayer.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/clusters.Plo@am__quote@
-@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/curvedelems_new.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/curvedelems.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/delaunay.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/geomsearch.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/global.Plo@am__quote@

libsrc/meshing/curvedelems.hpp

@@ -3,866 +3,207 @@
 
 /**************************************************************************/
 /* File:   curvedelems.hpp                                                */
-/* Author: Robert Gaisbauer                                               */
+/* Author: Robert Gaisbauer (first version)                               */
-/* Date:   27. Sep. 02 (second version: 30. Jan. 03)                      */
+/*         redesign by Joachim Schoeberl                                  */
+/* Date:   27. Sep. 02, Feb 2006                                          */
 /**************************************************************************/
 
-#include "bisect.hpp"
-#include <iostream>
-
-#define EPSILON 1e-20
 
 
 
-void ComputeGaussRule (int n, ARRAY<double> & xi, ARRAY<double> & wi);
+class Refinement;
 
 
-
-
-
-// ----------------------------------------------------------------------------
-//      CurvedElements
-// ----------------------------------------------------------------------------
-
 class CurvedElements
 {
   const Mesh & mesh;
-  const MeshTopology & top;
-
-  bool isHighOrder;
-  int nvisualsubsecs;
-  int nIntegrationPoints;
 
   ARRAY<int> edgeorder;
   ARRAY<int> faceorder;
 
-  /*
+  ARRAY<int> edgecoeffsindex;
+  ARRAY<int> facecoeffsindex;
 
   ARRAY< Vec<3> > edgecoeffs;
   ARRAY< Vec<3> > facecoeffs;
 
-  ARRAY<int> edgecoeffsindex;
+  ARRAY< double > edgeweight;  // for rational 2nd order splines
-  ARRAY<int> facecoeffsindex;
-
-  */
-
-  inline Vec<3> GetEdgeCoeff (int edgenr, int k);
-  inline Vec<3> GetFaceCoeff (int facenr, int k);
-
-  
-  void CalcSegmentTransformation (double xi, int segnr,
-				  Point<3> * x = NULL, Vec<3> * dxdxi = NULL);
-
-  void CalcSurfaceTransformation (Point<2> xi, int elnr,
-				  Point<3> * x = NULL, Mat<3,2> * dxdxi = NULL);
-
-  void CalcElementTransformation (Point<3> xi, int elnr,
-				  Point<3> * x = NULL, Mat<3,3> * dxdxi = NULL);
 
+  int order;
+  bool rational;
 
+  bool ishighorder;
 
 public:
-
-  Refinement * refinement;
-
-  ARRAY< Vec<3> > edgecoeffs;
-  ARRAY< Vec<3> > facecoeffs;
-
-  ARRAY<int> edgecoeffsindex;
-  ARRAY<int> facecoeffsindex;
-
-
-
-
-
   CurvedElements (const Mesh & amesh);
   ~CurvedElements();
 
-  bool IsHighOrder() const
+  // bool IsHighOrder() const { return order > 1; }
-  { return isHighOrder; };
+  bool IsHighOrder() const { return ishighorder; }
-  void SetHighOrder () { isHighOrder = 1; }
+
+  // void SetHighOrder (int aorder) { order=aorder; }
+  void SetIsHighOrder (bool ho) { ishighorder = ho; }
+  
+  void BuildCurvedElements(Refinement * ref, int aorder, bool arational = false);
+
+  int GetOrder () { return order; }
 
 
-  int GetNVisualSubsecs() const
+  bool IsSegmentCurved (SegmentIndex segnr) const;
-  { return nvisualsubsecs; };
+  bool IsSurfaceElementCurved (SurfaceElementIndex sei) const;
-
+  bool IsElementCurved (ElementIndex ei) const;
-  const class Mesh & GetMesh() const
-  { return mesh; };
-
-  void BuildCurvedElements(Refinement * ref, int polydeg, bool rational=false);
-
-  int GetEdgeOrder (int edgenr) const
-  { return edgeorder[edgenr]; };
-
-  int GetFaceOrder (int facenr) const
-  { return faceorder[facenr]; };
-
-  int IsEdgeCurved (int edgenr) const;
-
-  int IsFaceCurved (int facenr) const;
-
-  int IsSurfaceElementCurved (int elnr) const;
-
-  int IsElementCurved (int elnr) const;
 
 
-  void CalcSegmentTransformation (double xi, int segnr,
+  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
 				  Point<3> & x)
   { CalcSegmentTransformation (xi, segnr, &x, NULL); };
 
-  void CalcSegmentTransformation (double xi, int segnr,
+  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
 				  Vec<3> & dxdxi)
   { CalcSegmentTransformation (xi, segnr, NULL, &dxdxi); };
 
-  void CalcSegmentTransformation (double xi, int segnr,
+  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
 				  Point<3> & x, Vec<3> & dxdxi)
-  { CalcSegmentTransformation (xi, segnr, &x, &dxdxi); };
+  { CalcSegmentTransformation (xi, segnr, &x, &dxdxi, NULL); };
+
+  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
+				  Point<3> & x, Vec<3> & dxdxi, bool & curved)
+  { CalcSegmentTransformation (xi, segnr, &x, &dxdxi, &curved); };
 
 
-  void CalcSurfaceTransformation (const Point<2> & xi, int elnr,
+
+  void CalcSurfaceTransformation (const Point<2> & xi, SurfaceElementIndex elnr,
 				  Point<3> & x)
   { CalcSurfaceTransformation (xi, elnr, &x, NULL); };
 
-  void CalcSurfaceTransformation (const Point<2> & xi, int elnr,
+  void CalcSurfaceTransformation (const Point<2> & xi, SurfaceElementIndex elnr,
 				  Mat<3,2> & dxdxi)
   { CalcSurfaceTransformation (xi, elnr, NULL, &dxdxi); };
 
-  void CalcSurfaceTransformation (const Point<2> & xi, int elnr,
+  void CalcSurfaceTransformation (const Point<2> & xi, SurfaceElementIndex elnr,
 				  Point<3> & x, Mat<3,2> & dxdxi)
-  { CalcSurfaceTransformation (xi, elnr, &x, &dxdxi); };
+  { CalcSurfaceTransformation (xi, elnr, &x, &dxdxi, NULL); };
+
+  void CalcSurfaceTransformation (const Point<2> & xi, SurfaceElementIndex elnr,
+				  Point<3> & x, Mat<3,2> & dxdxi, bool & curved)
+  { CalcSurfaceTransformation (xi, elnr, &x, &dxdxi, &curved); };
 
 
-  void CalcElementTransformation (const Point<3> & xi, int elnr,
+
+
+
+  void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
 				  Point<3> & x)
   { CalcElementTransformation (xi, elnr, &x, NULL); };
 
-  void CalcElementTransformation (const Point<3> & xi, int elnr,
+  void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
 				  Mat<3,3> & dxdxi)
   { CalcElementTransformation (xi, elnr, NULL, &dxdxi); };
 
-  void CalcElementTransformation (const Point<3> & xi, int elnr,
+  void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
 				  Point<3> & x, Mat<3,3> & dxdxi)
-  { CalcElementTransformation (xi, elnr, &x, &dxdxi); };
+  { CalcElementTransformation (xi, elnr, &x, &dxdxi /* , NULL */ ); };
+
+  void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
+				  Point<3> & x, Mat<3,3> & dxdxi,
+                                  void * buffer, bool valid)
+  { CalcElementTransformation (xi, elnr, &x, &dxdxi, /* NULL, */ buffer, valid ); };
+
+  // void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
+  // 				  Point<3> & x, Mat<3,3> & dxdxi) // , bool & curved)
+  //   { CalcElementTransformation (xi, elnr, &x, &dxdxi /* , &curved * ); }
 
 
 
-
+  void CalcMultiPointSegmentTransformation (ARRAY<double> * xi, SegmentIndex segnr,
-  void CalcMultiPointSegmentTransformation (ARRAY<double> * xi, int segnr,
 					    ARRAY<Point<3> > * x,
 					    ARRAY<Vec<3> > * dxdxi);
 
-  void CalcMultiPointSurfaceTransformation (ARRAY< Point<2> > * xi, int elnr,
+  void CalcMultiPointSurfaceTransformation (ARRAY< Point<2> > * xi, SurfaceElementIndex elnr,
 					    ARRAY< Point<3> > * x,
 					    ARRAY< Mat<3,2> > * dxdxi);
 
-  void CalcMultiPointElementTransformation (ARRAY< Point<3> > * xi, int elnr,
+  void CalcMultiPointElementTransformation (ARRAY< Point<3> > * xi, ElementIndex elnr,
 					    ARRAY< Point<3> > * x,
 					    ARRAY< Mat<3,3> > * dxdxi);
 
-};
+  void CalcMultiPointElementTransformation (ElementIndex elnr, int n,
+                                            const double * xi, int sxi,
+                                            double * x, int sx,
+                                            double * dxdxi, int sdxdxi);
 
 
 
-// ----------------------------------------------------------------------------
-//      PolynomialBasis
-// ----------------------------------------------------------------------------
 
-class PolynomialBasis
+private:
-{
+  
-  int order;
+  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
-  int maxorder;
+				  Point<3> * x = NULL, Vec<3> * dxdxi = NULL, bool * curved = NULL);
-  ArrayMem<double,20> f;
-  ArrayMem<double,20> df;
-  ArrayMem<double,20> ddf;
 
-  ArrayMem<double,20> lp;
+  void CalcSurfaceTransformation (Point<2> xi, SurfaceElementIndex elnr,
-  ArrayMem<double,20> dlp;
+				  Point<3> * x = NULL, Mat<3,2> * dxdxi = NULL, bool * curved = NULL);
 
-  inline void CalcLegendrePolynomials (double x);
+  void CalcElementTransformation (Point<3> xi, ElementIndex elnr,
-  // P_i(x/t) t^i
+				  Point<3> * x = NULL, Mat<3,3> * dxdxi = NULL, // bool * curved = NULL,
-  inline void CalcScaledLegendrePolynomials (double x, double t);
+                                  void * buffer = NULL, bool valid = 0);
-  inline void CalcDLegendrePolynomials (double x);
 
-public:
 
-  PolynomialBasis ()
-  { maxorder = -1; };
 
-  ~PolynomialBasis ()
-  {};
 
-  void SetOrder (int aorder)
+
+
+  class SegmentInfo
   {
-    order = aorder;
+  public:
-    if (order > maxorder)
+    SegmentIndex elnr;
-      {
+    int order;
-	maxorder = order;
+    int nv;
-	f.SetSize(order-1);
+    int ndof;
-	df.SetSize(order-1);
+    int edgenr;
-	ddf.SetSize(order-1);
-	lp.SetSize(order+1);
-	dlp.SetSize(order);
-      };
   };
 
-  inline void CalcF (double x);
+  void CalcElementShapes (SegmentInfo &  elnr, double xi, Vector & shapes) const;
-  inline void CalcDf (double x);
+  void GetCoefficients (SegmentInfo & elnr, ARRAY<Vec<3> > & coefs) const;
-  inline void CalcDDf (double x);
+  void CalcElementDShapes (SegmentInfo & elnr, double xi, Vector & dshapes) const;
-
-  inline void CalcFDf (double x);
-
-  // compute F_i(x/t) t^i
-  inline void CalcFScaled (double x, double t);
-  static inline void CalcFScaled (int p, double x, double t, double * values);
-
-  double GetF (int p) { return f[p-2]; };
-  double GetDf (int p) { return df[p-2]; };
-  double GetDDf (int p) { return ddf[p-2]; };
-};
 
 
-
+  class ElementInfo
-// ----------------------------------------------------------------------------
+  {
-//      BaseFiniteElement
+  public:
-// ----------------------------------------------------------------------------
+    ElementIndex elnr;
-
+    int order;
-template <int DIM>
+    int nv;
-class BaseFiniteElement
+    int ndof;
-{
+    int nedges;
-protected:
+    int nfaces;
-
+    int edgenrs[12];
-  Point<DIM> xi;
+    int facenrs[6];
-  int elnr;
+    Mat<3> hdxdxi;
-  const CurvedElements & curv;
+    Vec<3> hcoefs[10]; // enough for second order tets
-  const Mesh & mesh;
-  const MeshTopology & top;
-
-public:
-
-  BaseFiniteElement(const CurvedElements & acurv)
-    : curv(acurv), mesh(curv.GetMesh()), top(mesh.GetTopology())
-  {};
-
-  virtual ~BaseFiniteElement()
-  {};
-
-  void SetElementNumber (int aelnr)
-  { elnr = aelnr; }; // 1-based arrays in netgen
-
-  virtual void SetReferencePoint (Point<DIM> axi)
-  { xi = axi; };
-};
-
-
-
-// ----------------------------------------------------------------------------
-//      BaseFiniteElement1D
-// ----------------------------------------------------------------------------
-
-class BaseFiniteElement1D : public BaseFiniteElement<1>
-{
-protected:
-  PolynomialBasis b;
-
-  int vertexnr[2];
-  int edgenr;
-  int edgeorient;
-  int edgeorder;
-
-  int maxedgeorder;
-
-  double vshape[2];
-  double vdshape[2];
-  ArrayMem<double,20> eshape;
-  ArrayMem<double,20> edshape;
-  ArrayMem<double,20> eddshape;
-
-public:
-
-  BaseFiniteElement1D (const CurvedElements & acurv) : BaseFiniteElement<1>(acurv)
-  { maxedgeorder = 1; };
-
-  virtual ~BaseFiniteElement1D()
-  {};
-
-  int GetVertexNr (int v)
-  { return vertexnr[v]; };
-
-  int GetEdgeNr ()
-  { return edgenr; };
-
-  int GetEdgeOrder ()
-  { return edgeorder; };
-
-  int GetEdgeOrientation ()
-  { return edgeorient; };
-
-  void CalcVertexShapes();
-  void CalcEdgeShapes();
-  void CalcEdgeLaplaceShapes();
-
-  double GetVertexShape (int v)
-  { return vshape[v]; };
-
-  double GetEdgeShape (int index)
-  { return eshape[index]; };
-
-  double GetVertexDShape (int v)
-  { return vdshape[v]; };
-
-  double GetEdgeDShape (int index)
-  { return edshape[index]; };
-
-  double GetEdgeLaplaceShape (int index)
-  { return eddshape[index]; };
-
-};
-
-
-
-
-// ----------------------------------------------------------------------------
-//      FESegm
-// ----------------------------------------------------------------------------
-
-class FESegm : public BaseFiniteElement1D
-{
-
-public:
-
-  FESegm(const CurvedElements & acurv) : BaseFiniteElement1D(acurv)
-  {};
-
-  virtual ~FESegm()
-  {};
-
-  void SetElementNumber (int aelnr)
-  { 
-    BaseFiniteElement<1> :: SetElementNumber (aelnr);
-    Segment s = mesh.LineSegment(elnr);
-    vertexnr[0] = s.p1;
-    vertexnr[1] = s.p2;
-    edgenr = top.GetSegmentEdge(elnr);
-    edgeorient = top.GetSegmentEdgeOrientation(elnr);
-    edgeorder = curv.GetEdgeOrder(edgenr-1); // 1-based arrays in netgen
-
-    if (edgeorder > maxedgeorder)
-      {
-	maxedgeorder = edgeorder;
-	eshape.SetSize(maxedgeorder-1);
-	edshape.SetSize(maxedgeorder-1);
-	eddshape.SetSize(maxedgeorder-1);
-      }
   };
 
+
+  void CalcElementShapes (ElementInfo & info, const Point<3> & xi, Vector & shapes) const;
+  void GetCoefficients (ElementInfo & info, Vec<3> * coefs) const;
+  void CalcElementDShapes (ElementInfo & info, const Point<3> & xi, MatrixFixWidth<3> & dshapes) const;
+
+  
+  class SurfaceElementInfo
+  {
+  public:
+    SurfaceElementIndex elnr;
+    int order;
+    int nv;
+    int ndof;
+    ArrayMem<int,4> edgenrs;
+    int facenr;
+  };
+
+  void CalcElementShapes (SurfaceElementInfo & elinfo, const Point<2> & xi, Vector & shapes) const;
+  void GetCoefficients (SurfaceElementInfo & elinfo, ARRAY<Vec<3> > & coefs) const;
+  void CalcElementDShapes (SurfaceElementInfo & elinfo, const Point<2> & xi, DenseMatrix & dshapes) const;
 };
 
 
 
-// ----------------------------------------------------------------------------
-//      FEEdge
-// ----------------------------------------------------------------------------
-
-class FEEdge : public BaseFiniteElement1D
-{
-
-public:
-
-  FEEdge(const CurvedElements & acurv) : BaseFiniteElement1D(acurv)
-  {};
-
-  virtual ~FEEdge()
-  {};
-
-  void SetElementNumber (int aelnr)
-  { 
-    BaseFiniteElement<1> :: SetElementNumber (aelnr);
-    top.GetEdgeVertices (elnr, vertexnr[0], vertexnr[1]);
-    edgenr = elnr;
-    edgeorient = 1;
-    edgeorder = curv.GetEdgeOrder(edgenr-1); // 1-based arrays in netgen
-
-    if (edgeorder > maxedgeorder)
-      {
-	maxedgeorder = edgeorder;
-	eshape.SetSize(maxedgeorder-1);
-	edshape.SetSize(maxedgeorder-1);
-	eddshape.SetSize(maxedgeorder-1);
-      }
-  };
-    
-};
-
-
-
-// ----------------------------------------------------------------------------
-//      BaseFiniteElement2D
-// ----------------------------------------------------------------------------
-
-class BaseFiniteElement2D : public BaseFiniteElement<2>
-{
-protected:
-
-  int nvertices;
-  int nedges;
-
-  int vertexnr[4];
-  int edgenr[4];
-  int edgeorient[4];
-  int edgeorder[4];
-  int facenr;
-  int faceorient;
-  int faceorder;
- 
-  int nfaceshapes;
-
-  int maxedgeorder;
-  int maxfaceorder;
-
-  PolynomialBasis b1, b2;
-
-  double vshape[4];
-  Vec<2> vdshape[4];
-  ArrayMem<double,80> eshape;
-  ArrayMem< Vec<2>,80> edshape;
-  ArrayMem<double,400> fshape;
-  ArrayMem<Vec<2>,400> fdshape;
-  ArrayMem<double,400> fddshape;
-
-  virtual void CalcNFaceShapes () = 0;
-
-public:
-
-  BaseFiniteElement2D (const CurvedElements & acurv) : BaseFiniteElement<2>(acurv)
-  { maxedgeorder = maxfaceorder = -1; };
-
-    virtual ~BaseFiniteElement2D()
-	{};
-
-  void SetElementNumber (int aelnr);
-
-  virtual void SetVertexSingularity (int v, int exponent) = 0;
-
-  int GetVertexNr (int v)
-  { return vertexnr[v]; };
-
-  int GetEdgeNr (int e)
-  { return edgenr[e]; };
-
-  int GetFaceNr ()
-  { return facenr; };
-
-  int GetEdgeOrder (int e)
-  { return edgeorder[e]; };
-
-  int GetFaceOrder ()
-  { return faceorder; }
-
-  int GetNVertices ()
-  { return nvertices; };
-
-  int GetNEdges ()
-  { return nedges; };
-
-  int GetNFaceShapes ()
-  { return nfaceshapes; };
-
-  int IsCurved ()
-  {
-    bool iscurved = 0;
-    int e;
-
-    for (e = 0; e < GetNEdges(); e++)
-      iscurved = iscurved || (GetEdgeOrder(e) > 1);
-
-    return iscurved || (GetFaceOrder() > 1);
-  }
-
-  virtual void CalcVertexShapes() = 0;
-  virtual void CalcEdgeShapes() = 0; 
-  virtual void CalcFaceShapes() = 0;
-
-  virtual void CalcFaceLaplaceShapes() = 0;
-
-  double GetVertexShape (int v)
-  { return vshape[v]; };
-
-  double GetEdgeShape (int index)
-  { return eshape[index]; };
-
-  double GetFaceShape (int index)
-  { return fshape[index]; };
-
-  Vec<2> GetVertexDShape (int v)
-  { return vdshape[v]; };
-
-  Vec<2> GetEdgeDShape (int index)
-  { return edshape[index]; };
-
-  Vec<2> GetFaceDShape (int index)
-  { return fdshape[index]; };
-
-  double GetFaceLaplaceShape (int index)
-  { return fddshape[index]; };
-};
-
-
-
-// ----------------------------------------------------------------------------
-//      FETrig
-// ----------------------------------------------------------------------------
-
-class FETrig : public BaseFiniteElement2D
-{
-  Point<3> lambda;
-  Mat<3,2> dlambda;
-
-  const ELEMENT_EDGE * eledge;
-  const ELEMENT_FACE * elface;
-
-  virtual void CalcNFaceShapes ()
-  { nfaceshapes = ((faceorder-1)*(faceorder-2))/2; };
-
-public:
-
-  FETrig (const CurvedElements & acurv) : BaseFiniteElement2D(acurv)
-  {
-    nvertices = 3;
-    nedges = 3;
-    eledge = MeshTopology :: GetEdges (TRIG);
-    elface = MeshTopology :: GetFaces (TRIG);
-  };
-
-    virtual ~FETrig()
-	{};
-
-  virtual void SetReferencePoint (Point<2> axi);
-
-  virtual void SetVertexSingularity (int v, int exponent);
-
-  virtual void CalcVertexShapes();
-  virtual void CalcEdgeShapes();
-  virtual void CalcFaceShapes();
-
-  virtual void CalcFaceLaplaceShapes();
-};
-
-
-
-// ----------------------------------------------------------------------------
-//      FEQuad
-// ----------------------------------------------------------------------------
-
-class FEQuad : public BaseFiniteElement2D
-{
-  const ELEMENT_FACE * elface;
-
-  virtual void CalcNFaceShapes ()
-  { nfaceshapes = (faceorder-1)*(faceorder-1); };
-
-public:
-
-  FEQuad (const CurvedElements & acurv) : BaseFiniteElement2D(acurv)
-  {
-    nvertices = 4;
-    nedges = 4;
-    elface = MeshTopology :: GetFaces (QUAD);
-  };
-
-    virtual ~FEQuad()
-	{};
-
-  virtual void SetVertexSingularity (int /* v */, int /* exponent */)
-	{};
-
-  virtual void CalcVertexShapes();
-  virtual void CalcEdgeShapes();
-  virtual void CalcFaceShapes();
-
-  virtual void CalcFaceLaplaceShapes();
-};
-
-
-
-
-// ----------------------------------------------------------------------------
-//      BaseFiniteElement3D
-// ----------------------------------------------------------------------------
-
-class BaseFiniteElement3D : public BaseFiniteElement<3>
-{
-protected:
-
-  int nvertices;
-  int nedges;
-  int nfaces;
-
-  int vertexnr[8];
-  int edgenr[12];
-  int edgeorient[12];
-  int edgeorder[12];
-  int facenr[6];
-  int faceorient[6];
-  int faceorder[6];
-  int surfacenr[6];
-  // int surfaceorient[6];
-
-  int nfaceshapes[6];
-
-  int maxedgeorder;
-  int maxfaceorder;
-
-  PolynomialBasis b1, b2;
-
-  double vshape[8];
-  Vec<3> vdshape[8];
-  ArrayMem<double,120> eshape;
-  ArrayMem<Vec<3>,120> edshape;
-  ArrayMem<double,300> fshape;
-  ArrayMem<Vec<3>,300> fdshape;
-
-  virtual void CalcNFaceShapes () = 0;
-
-public:
-
-  int locmaxedgeorder;
-  int locmaxfaceorder;
-
-  BaseFiniteElement3D (const CurvedElements & acurv) : BaseFiniteElement<3>(acurv)
-  { maxedgeorder = maxfaceorder = -1; };
-
-  void SetElementNumber (int aelnr);
-
-  int GetVertexNr (int v)
-  { return vertexnr[v]; };
-
-  int GetEdgeNr (int e)
-  { return edgenr[e]; };
-
-  int GetFaceNr (int f)
-  { return facenr[f]; };
-
-  int GetNFaceShapes (int f)
-  { return nfaceshapes[f]; };
-
-  int GetEdgeOrder (int e)
-  { return edgeorder[e]; };
-
-  int GetFaceOrder (int f)
-  { return faceorder[f]; };
-
-  int GetNVertices ()
-  { return nvertices; };
-
-  int GetNEdges ()
-  { return nedges; };
-
-  int GetNFaces ()
-  { return nfaces; };
-
-  int IsCurved ()
-  {
-    bool iscurved = 0;
-    int e, f;
-
-    for (e = 0; e < GetNEdges(); e++)
-      iscurved = iscurved || (GetEdgeOrder(e) > 1);
-
-    for (f = 0; f < GetNFaces(); f++)
-      iscurved = iscurved || (GetFaceOrder(f) > 1);
-
-    return iscurved;
-  }
-
-  virtual void CalcVertexShapes() = 0;
-  virtual void CalcVertexShapesOnly()
-  { CalcVertexShapes(); }
-
-  virtual void CalcEdgeShapes() = 0;
-  virtual void CalcEdgeShapesOnly() 
-    { CalcEdgeShapes(); }
-
-  virtual void CalcFaceShapes() = 0;
-
-  double GetVertexShape (int v)
-  { return vshape[v]; };
-
-  double GetEdgeShape (int index)
-  { return eshape[index]; };
-
-  double GetFaceShape (int index)
-  { return fshape[index]; };
-
-  Vec<3> GetVertexDShape (int v)
-  { return vdshape[v]; };
-
-  Vec<3> GetEdgeDShape (int index)
-  { return edshape[index]; };
-
-  Vec<3> GetFaceDShape (int index)
-  { return fdshape[index]; };
-};
-
-
-
-// ----------------------------------------------------------------------------
-//      FETet
-// ----------------------------------------------------------------------------
-
-class FETet : public BaseFiniteElement3D
-{
-  Point<4> lambda;
-  Mat<4,3> dlambda;
-
-  const ELEMENT_EDGE * eledge;
-  const ELEMENT_FACE * elface;
-
-  virtual void CalcNFaceShapes ()
-  {
-    for (int f = 0; f < nfaces; f++)
-      nfaceshapes[f] = ((faceorder[f]-1)*(faceorder[f]-2))/2;
-  };
-
-public:
-
-  FETet (const CurvedElements & acurv) : BaseFiniteElement3D(acurv)
-  {
-    nvertices = 4;
-    nedges = 6;
-    nfaces = 4;
-    eledge = MeshTopology :: GetEdges (TET);
-    elface = MeshTopology :: GetFaces (TET);
-  };
-
-  void SetReferencePoint (Point<3> axi);
-
-  virtual void CalcVertexShapes();
-  virtual void CalcVertexShapesOnly();
-  virtual void CalcEdgeShapes();
-  virtual void CalcEdgeShapesOnly();
-  virtual void CalcFaceShapes();
-};
-
-
-
-// ----------------------------------------------------------------------------
-//      FEPrism
-// ----------------------------------------------------------------------------
-
-class FEPrism : public BaseFiniteElement3D
-{
-  Point<4> lambda;   // mixed barycentric coordinates
-  Mat<4,3> dlambda;
-
-  const ELEMENT_EDGE * eledge;
-  const ELEMENT_FACE * elface;
-
-  virtual void CalcNFaceShapes ()
-  {
-    int f;
-    for (f = 0; f < 2; f++)
-      nfaceshapes[f] = ((faceorder[f]-1)*(faceorder[f]-2))/2;
-    for (f = 2; f < nfaces; f++)
-      nfaceshapes[f] = (faceorder[f]-1)*(faceorder[f]-1);
-  };
-
-public:
-
-  FEPrism (const CurvedElements & acurv) : BaseFiniteElement3D(acurv)
-  {
-    nvertices = 6;
-    nedges = 9;
-    nfaces = 5;
-    eledge = MeshTopology :: GetEdges (PRISM);
-    elface = MeshTopology :: GetFaces (PRISM);
-  };
-
-  void SetReferencePoint (Point<3> axi);
-
-  virtual void CalcVertexShapes();
-  virtual void CalcEdgeShapes();
-  virtual void CalcFaceShapes();
-};
-
-
-
-
-// ----------------------------------------------------------------------------
-//      FEPyramid
-// ----------------------------------------------------------------------------
-
-class FEPyramid : public BaseFiniteElement3D
-{
-
-  const ELEMENT_EDGE * eledge;
-  const ELEMENT_FACE * elface;
-
-  virtual void CalcNFaceShapes ()
-  {
-    int f;
-    for (f = 0; f < 4; f++)
-      nfaceshapes[f] = ((faceorder[f]-1)*(faceorder[f]-2))/2;
-    for (f = 4; f < nfaces; f++)
-      nfaceshapes[f] = (faceorder[f]-1)*(faceorder[f]-1);
-  };
-
-public:
-
-  FEPyramid (const CurvedElements & acurv) : BaseFiniteElement3D(acurv)
-  {
-    nvertices = 5;
-    nedges = 8;
-    nfaces = 5;
-    eledge = MeshTopology :: GetEdges (PYRAMID);
-    elface = MeshTopology :: GetFaces (PYRAMID);
-  };
-
-  void SetReferencePoint (Point<3> axi);
-
-  virtual void CalcVertexShapes();
-  virtual void CalcEdgeShapes();
-  virtual void CalcFaceShapes();
-};
-
-
-
-
-// ----------------------------------------------------------------------------
-//      FEHex
-// ----------------------------------------------------------------------------
-
-class FEHex : public BaseFiniteElement3D
-{
-
-  const ELEMENT_EDGE * eledge;
-  const ELEMENT_FACE * elface;
-
-  virtual void CalcNFaceShapes ()
-  {
-    int f;
-    for (f = 0; f < 6; f++)
-      nfaceshapes[f] = (faceorder[f]-1)*(faceorder[f]-1);
-  };
-
-public:
-
-  FEHex (const CurvedElements & acurv) : BaseFiniteElement3D(acurv)
-  {
-    nvertices = 8;
-    nedges = 12;
-    nfaces = 6;
-    eledge = MeshTopology :: GetEdges (HEX);
-    elface = MeshTopology :: GetFaces (HEX);
-  };
-
-  void SetReferencePoint (Point<3> axi);
-
-  virtual void CalcVertexShapes();
-  virtual void CalcEdgeShapes();
-  virtual void CalcFaceShapes();
-};
-
-
-
-
 #endif

libsrc/meshing/curvedelems_new.hpp

@@ -1,209 +0,0 @@
-#ifndef CURVEDELEMS_NEWH
-#define CURVEDELEMS_NEWH
-
-/**************************************************************************/
-/* File:   curvedelems.hpp                                                */
-/* Author: Robert Gaisbauer (first version)                               */
-/*         redesign by Joachim Schoeberl                                  */
-/* Date:   27. Sep. 02, Feb 2006                                          */
-/**************************************************************************/
-
-
-
-
-class Refinement;
-
-
-class CurvedElements
-{
-  const Mesh & mesh;
-
-  ARRAY<int> edgeorder;
-  ARRAY<int> faceorder;
-
-  ARRAY<int> edgecoeffsindex;
-  ARRAY<int> facecoeffsindex;
-
-  ARRAY< Vec<3> > edgecoeffs;
-  ARRAY< Vec<3> > facecoeffs;
-
-  ARRAY< double > edgeweight;  // for rational 2nd order splines
-
-  int order;
-  bool rational;
-
-  bool ishighorder;
-
-public:
-  CurvedElements (const Mesh & amesh);
-  ~CurvedElements();
-
-  // bool IsHighOrder() const { return order > 1; }
-  bool IsHighOrder() const { return ishighorder; }
-
-  // void SetHighOrder (int aorder) { order=aorder; }
-  void SetIsHighOrder (bool ho) { ishighorder = ho; }
-  
-  void BuildCurvedElements(Refinement * ref, int aorder, bool arational = false);
-
-  int GetOrder () { return order; }
-
-
-  bool IsSegmentCurved (SegmentIndex segnr) const;
-  bool IsSurfaceElementCurved (SurfaceElementIndex sei) const;
-  bool IsElementCurved (ElementIndex ei) const;
-
-
-  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
-				  Point<3> & x)
-  { CalcSegmentTransformation (xi, segnr, &x, NULL); };
-
-  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
-				  Vec<3> & dxdxi)
-  { CalcSegmentTransformation (xi, segnr, NULL, &dxdxi); };
-
-  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
-				  Point<3> & x, Vec<3> & dxdxi)
-  { CalcSegmentTransformation (xi, segnr, &x, &dxdxi, NULL); };
-
-  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
-				  Point<3> & x, Vec<3> & dxdxi, bool & curved)
-  { CalcSegmentTransformation (xi, segnr, &x, &dxdxi, &curved); };
-
-
-
-  void CalcSurfaceTransformation (const Point<2> & xi, SurfaceElementIndex elnr,
-				  Point<3> & x)
-  { CalcSurfaceTransformation (xi, elnr, &x, NULL); };
-
-  void CalcSurfaceTransformation (const Point<2> & xi, SurfaceElementIndex elnr,
-				  Mat<3,2> & dxdxi)
-  { CalcSurfaceTransformation (xi, elnr, NULL, &dxdxi); };
-
-  void CalcSurfaceTransformation (const Point<2> & xi, SurfaceElementIndex elnr,
-				  Point<3> & x, Mat<3,2> & dxdxi)
-  { CalcSurfaceTransformation (xi, elnr, &x, &dxdxi, NULL); };
-
-  void CalcSurfaceTransformation (const Point<2> & xi, SurfaceElementIndex elnr,
-				  Point<3> & x, Mat<3,2> & dxdxi, bool & curved)
-  { CalcSurfaceTransformation (xi, elnr, &x, &dxdxi, &curved); };
-
-
-
-
-
-  void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
-				  Point<3> & x)
-  { CalcElementTransformation (xi, elnr, &x, NULL); };
-
-  void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
-				  Mat<3,3> & dxdxi)
-  { CalcElementTransformation (xi, elnr, NULL, &dxdxi); };
-
-  void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
-				  Point<3> & x, Mat<3,3> & dxdxi)
-  { CalcElementTransformation (xi, elnr, &x, &dxdxi /* , NULL */ ); };
-
-  void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
-				  Point<3> & x, Mat<3,3> & dxdxi,
-                                  void * buffer, bool valid)
-  { CalcElementTransformation (xi, elnr, &x, &dxdxi, /* NULL, */ buffer, valid ); };
-
-  // void CalcElementTransformation (const Point<3> & xi, ElementIndex elnr,
-  // 				  Point<3> & x, Mat<3,3> & dxdxi) // , bool & curved)
-  //   { CalcElementTransformation (xi, elnr, &x, &dxdxi /* , &curved * ); }
-
-
-
-  void CalcMultiPointSegmentTransformation (ARRAY<double> * xi, SegmentIndex segnr,
-					    ARRAY<Point<3> > * x,
-					    ARRAY<Vec<3> > * dxdxi);
-
-  void CalcMultiPointSurfaceTransformation (ARRAY< Point<2> > * xi, SurfaceElementIndex elnr,
-					    ARRAY< Point<3> > * x,
-					    ARRAY< Mat<3,2> > * dxdxi);
-
-  void CalcMultiPointElementTransformation (ARRAY< Point<3> > * xi, ElementIndex elnr,
-					    ARRAY< Point<3> > * x,
-					    ARRAY< Mat<3,3> > * dxdxi);
-
-  void CalcMultiPointElementTransformation (ElementIndex elnr, int n,
-                                            const double * xi, int sxi,
-                                            double * x, int sx,
-                                            double * dxdxi, int sdxdxi);
-
-
-
-
-private:
-  
-  void CalcSegmentTransformation (double xi, SegmentIndex segnr,
-				  Point<3> * x = NULL, Vec<3> * dxdxi = NULL, bool * curved = NULL);
-
-  void CalcSurfaceTransformation (Point<2> xi, SurfaceElementIndex elnr,
-				  Point<3> * x = NULL, Mat<3,2> * dxdxi = NULL, bool * curved = NULL);
-
-  void CalcElementTransformation (Point<3> xi, ElementIndex elnr,
-				  Point<3> * x = NULL, Mat<3,3> * dxdxi = NULL, // bool * curved = NULL,
-                                  void * buffer = NULL, bool valid = 0);
-
-
-
-
-
-
-  class SegmentInfo
-  {
-  public:
-    SegmentIndex elnr;
-    int order;
-    int nv;
-    int ndof;
-    int edgenr;
-  };
-
-  void CalcElementShapes (SegmentInfo &  elnr, double xi, Vector & shapes) const;
-  void GetCoefficients (SegmentInfo & elnr, ARRAY<Vec<3> > & coefs) const;
-  void CalcElementDShapes (SegmentInfo & elnr, double xi, Vector & dshapes) const;
-
-
-  class ElementInfo
-  {
-  public:
-    ElementIndex elnr;
-    int order;
-    int nv;
-    int ndof;
-    int nedges;
-    int nfaces;
-    int edgenrs[12];
-    int facenrs[6];
-    Mat<3> hdxdxi;
-    Vec<3> hcoefs[10]; // enough for second order tets
-  };
-
-
-  void CalcElementShapes (ElementInfo & info, const Point<3> & xi, Vector & shapes) const;
-  void GetCoefficients (ElementInfo & info, Vec<3> * coefs) const;
-  void CalcElementDShapes (ElementInfo & info, const Point<3> & xi, MatrixFixWidth<3> & dshapes) const;
-
-  
-  class SurfaceElementInfo
-  {
-  public:
-    SurfaceElementIndex elnr;
-    int order;
-    int nv;
-    int ndof;
-    ArrayMem<int,4> edgenrs;
-    int facenr;
-  };
-
-  void CalcElementShapes (SurfaceElementInfo & elinfo, const Point<2> & xi, Vector & shapes) const;
-  void GetCoefficients (SurfaceElementInfo & elinfo, ARRAY<Vec<3> > & coefs) const;
-  void CalcElementDShapes (SurfaceElementInfo & elinfo, const Point<2> & xi, DenseMatrix & dshapes) const;
-};
-
-
-
-#endif

libsrc/meshing/global.cpp

@@ -7,7 +7,11 @@ namespace netgen
   // stringstream emptystr;
   // ostream * testout = &emptystr;
   // testout -> clear(ios::failbit);
-  ostream * testout = &cout;
+
+  // ostream * testout = &cout;
+  ostream * testout = new ostream(0);
+
+  // NetgenOutStream * testout = new NetgenOutStream;
 
   ostream * mycout = &cout;
   ostream * myerr = &cerr;

libsrc/meshing/global.hpp

@@ -19,9 +19,6 @@ extern void ResetTime ();
 ///
 extern int testmode;
 
-// extern ostream * testout;
-// extern AutoPtr<ostream> testout;
-
 /// calling parameters
 extern Flags parameters;
 

libsrc/meshing/meshclass.cpp

@@ -101,6 +101,9 @@ namespace netgen
 
   void Mesh :: DeleteMesh()
   {
+    NgLock lock(mutex);
+    lock.Lock();
+
     points.SetSize(0);
     segments.SetSize(0);
     surfelements.SetSize(0);
@@ -131,6 +134,7 @@ namespace netgen
     paralleltop = new ParallelMeshTopology (*this);
 #endif
 
+    lock.UnLock();
 
     timestamp = NextTimeStamp();
   }

libsrc/meshing/meshing.hpp

@@ -3,9 +3,6 @@
 
 
 
-#define CURVEDELEMS_NEW
-
-
 #include "../include/myadt.hpp"
 #include "../include/gprim.hpp"
 #include "../include/linalg.hpp"
@@ -56,12 +53,7 @@ namespace netgen
 #include "findip2.hpp"
 
 #include "topology.hpp"
-
-#ifdef CURVEDELEMS_NEW
-#include "curvedelems_new.hpp"
-#else
 #include "curvedelems.hpp"
-#endif
 #include "clusters.hpp"
 
 #ifdef _INCLUDE_MORE

libsrc/meshing/meshtool.cpp

@@ -153,7 +153,8 @@ namespace netgen
 
     (*testout) << endl;
     (*testout) << "Elements in qualityclasses:" << endl;
-    (*testout).precision(2);
+    // (*testout).precision(2);
+    (*testout) << setprecision(2);
     for (i = 1; i <= ncl; i++)
       {
 	(*testout) << setw(4) << double (i-1)/ncl << " - "
@@ -607,7 +608,7 @@ namespace netgen
     (*testout) << endl;
 
     (*testout) << "Volume elements in qualityclasses:" << endl;
-    (*testout).precision(2);
+    (*testout) << setprecision(2);
     for (i = 1; i <= ncl; i++)
       {
 	(*testout) << setw(4) << double (i-1)/ncl << " - "

libsrc/meshing/smoothing3.cpp

@@ -1266,7 +1266,7 @@ void Mesh :: ImproveMesh (const CSGeometry & geometry, OPTIMIZEGOAL goal)
 
   int uselocalh = mparam.uselocalh;
 
-  (*testout).precision(8);
+  (*testout) << setprecision(8);
   (*testout) << "Improve Mesh" << "\n";
   PrintMessage (3, "ImproveMesh");
   //  (*mycout) << "Vol = " << CalcVolume (points, volelements) << endl;
@@ -1470,7 +1470,7 @@ void Mesh :: ImproveMesh (OPTIMIZEGOAL goal)
   
   Vector x(3);
   
-  (*testout).precision(8);
+  (*testout) << setprecision(8);
   
   //int uselocalh = mparam.uselocalh;
 
@@ -1606,7 +1606,7 @@ void Mesh :: ImproveMeshJacobian (OPTIMIZEGOAL goal, const BitArray * usepoint)
   
   Vector x(3);
   
-  (*testout).precision(8);
+  (*testout) << setprecision(8);
   
   JacobianPointFunction pf(points, volelements);
   

libsrc/visualization/vsmesh.cpp

@@ -330,7 +330,13 @@ namespace netgen
 	VisualScene::BuildScene (zoomall);
 	return;
       }
-      
+
+    if (!lock)
+      {
+	lock = new NgLock (mesh->Mutex());
+	lock -> Lock();
+      }
+
     int i, j;
 	
 
@@ -414,7 +420,6 @@ namespace netgen
 	identifiedlist = 0;
       }
 
-
     pointnumberlist = glGenLists (1);
     glNewList (pointnumberlist, GL_COMPILE);
 
@@ -604,13 +609,6 @@ namespace netgen
 
 
 
-
-
-
-
-
-
-
     badellist = glGenLists (1);
     glNewList (badellist, GL_COMPILE);
 
@@ -845,9 +843,6 @@ namespace netgen
     glEndList ();
   
 
-
-  
-
     if (1)
       {
       
@@ -858,33 +853,45 @@ namespace netgen
   
 	glLineWidth (3);
       
-
 	//  for (i = 1; i <= mesh->GetNSeg(); i++)
-	INDEX_2_HASHTABLE<int> & idpts = 
-	  mesh->GetIdentifications().GetIdentifiedPoints();
-	if (&idpts)
-	  for (i = 1; i <= idpts.GetNBags(); i++)
-	    for (j = 1; j <= idpts.GetBagSize(i); j++)
-	      {
-		INDEX_2 pts;
-		int val;
 
-		idpts.GetData (i, j, pts, val);
+        if (& mesh -> GetIdentifications() )
-		const Point3d & p1 = mesh->Point(pts.I1());
+          {
-		const Point3d & p2 = mesh->Point(pts.I2());
+            INDEX_2_HASHTABLE<int> & idpts = 
-      
+              mesh->GetIdentifications().GetIdentifiedPoints();
-		glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, 
+            if (&idpts)
-			      identifiedcol);
+              {
-
+                for (i = 1; i <= idpts.GetNBags(); i++)
-		glBegin (GL_LINES);
+                  for (j = 1; j <= idpts.GetBagSize(i); j++)
-		glVertex3f (p1.X(), p1.Y(), p1.Z());
+                    {
-		glVertex3f (p2.X(), p2.Y(), p2.Z());
+                      INDEX_2 pts;
-		glEnd(); 
+                      int val;
-	      }
+                      
+                      idpts.GetData (i, j, pts, val);
+                      const Point3d & p1 = mesh->Point(pts.I1());
+                      const Point3d & p2 = mesh->Point(pts.I2());
+                      
+                      glMaterialfv (GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, 
+                                    identifiedcol);
+                      
+                      glBegin (GL_LINES);
+                      glVertex3f (p1.X(), p1.Y(), p1.Z());
+                      glVertex3f (p2.X(), p2.Y(), p2.Z());
+                      glEnd(); 
+                    }
+              }
+          }
 
 	glEndList ();
       }
 
+    if (lock)
+      {
+	lock -> UnLock();
+	delete lock;
+	lock = NULL;
+      }
+
     vstimestamp = meshtimestamp;
   }
 
@@ -936,7 +943,6 @@ namespace netgen
     // clock_t starttime, endtime;
     // starttime = clock();
 
-    
     if (!lock)
       {
 	lock = new NgLock (mesh->Mutex());

ng/drawing.tcl

@@ -4,6 +4,7 @@
 set oldmousex 0
 set oldmousey 0
 #
+
 # if { 1 } {
 if {[catch {togl .ndraw -width 400 -height 300  -rgba true -double true -depth true -privatecmap false -stereo false -indirect false }] } {    
     puts "no OpenGL" 

ng/ng.tcl

@@ -89,7 +89,7 @@ catch {
 
 
 
-catch { source ${ngdir}/ngsolve/ngsolve.tcl }
+# catch { source ${ngdir}/ngsolve/ngsolve.tcl }
 # catch { [source ${ngdir}/ngsolve/preproc.tcl] }
 # catch { [source ${ngdir}/ngsolve/pdecreator.tcl] }
 

ng/ngappinit.cpp

@@ -183,8 +183,8 @@ int main(int argc, char ** argv)
 
   if ( netgen::id == 0 )
     {
-      // if (parameters.StringFlagDefined ("testout"))
+      if (parameters.StringFlagDefined ("testout"))      
-      netgen::testout = new ofstream (parameters.GetStringFlag ("testout", "test.out"));
+        netgen::testout = new ofstream (parameters.GetStringFlag ("testout", "test.out"));
 
 
 #ifdef SOCKETS
@@ -330,7 +330,7 @@ int main(int argc, char ** argv)
       MPI_Barrier(MPI_COMM_WORLD);
       MPI_Finalize();
 #endif
-
+      
       Tcl_Exit(0);
     }
 #ifdef PARALLEL

ng/ngpkg.cpp

@@ -4596,7 +4596,9 @@ namespace netgen
     for ( int dest = 1; dest < ntasks; dest++)
       MyMPI_Send ( "end", dest );
 #endif
-    delete testout;
+
+    if (testout != &cout)
+      delete testout;
 
     return TCL_OK;
   }

nglib/nglib.cpp

@@ -46,6 +46,7 @@ void Ng_Init ()
 {
   mycout = &cout;
   myerr = &cerr;
+  // netgen::testout->SetOutStream (new ofstream ("test.out"));
   testout = new ofstream ("test.out");
 }