I/O for the amiramesh format; no functionality so far

[[Imported from SVN: r357]]

lib/io/amirameshreader.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef __DUNE_AMIRAMESHREADER_HH__
+#define __DUNE_AMIRAMESHREADER_HH__
+
+#include <string>
+#include "../../common/array.hh"
+
+namespace Dune {
+
+  /** @ingroup IO
+   * \brief Provides file reading facilities in the AmiraMesh format.
+   *
+   */
+  template<class GRID, class T>
+  class AmiraMeshReader {
+
+  public:
+
+    /** \brief The method that does the writing.
+     *
+     * @param grid The grid objects that is to be written
+     * @param sol  Data that should be written along with the grid
+     * @param filename The filename
+     */
+    static void read(GRID& grid,
+                     const std::string& filename);
+
+
+    AmiraMeshReader() {}
+
+  };
+
+}
+
+// Default implementation
+template<class GRID, class T>
+void Dune::AmiraMeshReader<GRID, T>::read(GRID& grid,
+                                          const std::string& filename)
+{
+  printf("No AmiraMesh reading has been implemented for this grid type!\n");
+}
+
+
+// the amiramesh reader for UGGrid
+#ifdef HAVE_UG
+#include "amuggridreader.cc"
+#endif
+
+#endif

lib/io/amirameshwriter.hh

@@ -50,6 +50,15 @@ void Dune::AmiraMeshWriter<GRID, T>::write(const GRID& grid,
   printf("No AmiraMesh writing has been implemented for this grid type!\n");
 }
 
+
+// the amiramesh writer for SimpleGrid
+#ifndef __GNUC__
 #include "amsimplegridwriter.cc"
+#endif
+
+// the amiramesh writer for UGGrid
+#ifdef HAVE_UG
+#include "amuggridwriter.cc"
+#endif
 
 #endif

lib/io/amuggridreader.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+// ///////////////////////////////////////////////
+// Specialization of the AmiraMesh reader for UGGrid<3,3>
+// ///////////////////////////////////////////////
+
+#include "../../grid/uggrid.hh"
+
+#include <amiramesh/AmiraMesh.h>
+#include "../../../AmiraLibs/include/AmiraParamAccess.h"
+
+namespace Dune {
+
+  template<>
+  class AmiraMeshReader<UGGrid<3,3>, double> {
+
+  public:
+
+    static void read(UGGrid<3,3>& grid,
+                     const std::string& filename);
+
+    static int CreateDomain(UGGrid<3,3>& grid,
+                            const std::string& domainName,
+                            const std::string& filename);
+
+    AmiraMeshReader() {}
+
+  };
+
+}
+
+
+// //////////////////////////////////////////////////
+// //////////////////////////////////////////////////
+static int SegmentDescriptionByAmira(void *data, double *param, double *result)
+{
+
+  int triNum = * (int*) data;
+
+  //printf("tri = %d, l0 = %5.5f, l1 = %5.5f\n", triNum, param[0], param[1]);
+
+  double barCoords[2];
+  double A[4] = {-1, 1, 0, -1};
+  double b[2] = {1, 0};
+
+  // barCoords = A*param + b;
+  barCoords[0] = A[0]*param[0] + A[2]*param[1];
+  barCoords[1] = A[1]*param[0] + A[3]*param[1];
+
+  barCoords[0] += b[0];
+  barCoords[1] += b[1];
+
+  /*
+     An dieser Stelle wird vorausgesetzt, dass das bereits gesetzte Gebiet
+     nicht mehr geaendert wird
+   */
+  //printf("barCoords:  (%g %g)\n", barCoords[0], barCoords[1]);
+
+  const float eps = 1e-6;
+  if (barCoords[0]<-eps || barCoords[1]<-eps || (1-barCoords[0]-barCoords[1])<-eps) {
+    printf("Illegal barycentric coordinate\n");
+    assert(false);
+  }
+
+  AmiraCallPositionParametrization(triNum, barCoords, result);
+
+  return(0);
+}
+
+
+int Dune::AmiraMeshReader<Dune::UGGrid<3,3>, double>::CreateDomain(UGGrid<3,3>& grid,
+                                                                   const std::string& domainName,
+                                                                   const std::string& filename)
+{
+  const int DIM = 3;
+  const int MAX_CORNERS_OF_LINEAR_PATCH = DIM;
+  const int CORNERS_OF_BND_SEG = 4;
+
+
+  int point[CORNERS_OF_BND_SEG] = {-1, -1, -1, -1};
+  double radius, MidPoint[3], alpha[2], beta[2], pos[3];
+
+  /* Alle Aufrufe an Deine Bibliothek beginnen mit "Amira" */
+
+  /* Load data */
+  if(AmiraLoadMesh(domainName.c_str(), filename.c_str()) != AMIRA_OK)
+  {
+    UG3d::PrintErrorMessage('E', "CreateAmiraParametrization", "Domain not found");
+    return(1);
+  }
+
+  if(AmiraStartEditingDomain(domainName.c_str()) != AMIRA_OK)
+  {
+    UG3d::PrintErrorMessage('E', "AmiraStartEditing", "Domain not found");
+    return(1);
+  }
+
+
+  /* Alle weiteren Anfragen an die Bibliothek beziehen sich jetzt auf das eben
+     geladen Gebiet. Maessig elegant, sollte aber gehen */
+
+  int noOfSegments = AmiraGetNoOfSegments();
+  if(noOfSegments <= 0)
+  {
+    UG3d::PrintErrorMessage('E', "CreateAmiraDomain", "no segments found");
+    return(1);
+  }
+
+  int noOfNodes = AmiraGetNoOfNodes();
+  if(noOfNodes <= 0)
+  {
+    UG3d::PrintErrorMessage('E', "CreateAmiraDomain", "no nodes found");
+    return(1);
+  }
+
+
+  /* Das naechste ist neu. Wir brauchen eine Kugel, die das Gebiet komplett  enthaelt.
+     Diese Information wird fuer die UG-Graphik benoetigt, die Kugel kann also durchaus
+     zu gross oder zu klein sein */
+
+  AmiraGetMidpointAndRadius(MidPoint, &radius);
+
+  //     printf("Amira radius = %5.2f, Midpoint = (%5.2f, %5.2f, %5.2f)\n",
+  //            radius, MidPoint[0], MidPoint[1], MidPoint[2]);
+
+
+  /* Jetzt geht es an's eingemachte. Zuerst wird ein neues gebiet konstruiert und
+     in der internen UG Datenstruktur eingetragen */
+
+  bool isConvex = false;
+
+  UG3d::domain* newDomain = UG3d::CreateDomain(const_cast<char*>(domainName.c_str()),
+                                               MidPoint, radius,
+                                               noOfSegments, noOfNodes,
+                                               isConvex );
+
+
+  if (!newDomain)
+    return(1);
+
+
+  /* Alle weiteren Aufrufe von 'CreateBoundarySegment' beziehen sich jetzt auf das eben
+     erzeugte Gebiet */
+
+
+  /*
+     alpha und beta spannen den  Parameterbereich auf (s. beiliegendes .ps-file)
+     wenn ich das richtig verstehe, sollten die konstant sein fuer alle Segmente
+   */
+
+  alpha[0] = 0.0;
+  alpha[1] = 0.0;
+  beta[0]  = 1.0;
+  beta[1]  = 1.0;
+
+  /*
+     Liste der Nummern der Randsegmente herstellen, wird als user data an das jeweilige
+     Segment weitergereicht
+   */
+
+  int* segmentIndexList = NULL;
+  segmentIndexList = (int*) ::malloc(noOfSegments*sizeof(int));
+
+  if(segmentIndexList == NULL)
+    return(1);
+
+  for(int i = 0; i < noOfSegments; i++) {
+
+    //std::string segmentName;
+    char segmentName[200];
+    int left, right, j;
+    double coordinates[MAX_CORNERS_OF_LINEAR_PATCH][DIM];
+
+    AmiraGetNodeNumbersOfSegment(point, i);
+
+    /* Es werden die Ecknummern eines Randsegmentes zurueckgegeben
+       point[0] = 0; point[1]=1; point[2]=2; point[3]=3; */
+
+    //segmentName = "AmiraSegment";
+    if(sprintf(segmentName, "AmiraSegment %d", i) < 0)
+      return(1);
+
+    /* left = innerRegion, right = outerRegion */
+    AmiraGetLeftAndRightSideOfSegment(&left, &right, i);
+
+    //             printf("id = %d l = %d, r = %d\n", i, left, right);
+    //             printf("id = %d nd0 = %d, nd1 = %d, nd2 = %d\n", i, point[0], point[1], point[2]);
+
+    segmentIndexList[i] = i;
+
+    /* 10 / 01 / 00 */
+    alpha[0] = 1.0; alpha[1] = 0.0;
+    AmiraCallPositionParametrization(i, alpha, pos);
+    for(j=0; j<DIM; j++)
+      coordinates[0][j] = pos[j];
+
+    alpha[0] = 0.0; alpha[1] = 1.0;
+    AmiraCallPositionParametrization(i, alpha, pos);
+    for(j=0; j<DIM; j++)
+      coordinates[1][j] = pos[j];
+
+    alpha[0] = 0.0; alpha[1] = 0.0;
+    AmiraCallPositionParametrization(i, alpha, pos);
+    for(j=0; j<DIM; j++)
+      coordinates[2][j] = pos[j];
+
+
+    /* map Amira Material ID's to UG material ID's */
+
+    left++;
+    right++;
+
+    alpha[0] = alpha[1] = 0;
+    beta[0]  = beta[1]  = 1;
+
+    if (UG3d::CreateBoundarySegment(segmentName,
+                                    left,           /*id of left subdomain */
+                                    right,          /*id of right subdomain*/
+                                    i,              /*id of segment*/
+                                    3,              // it' a triangle
+                                    1,              // resolution, whatever that is
+                                    point,
+                                    alpha, beta,
+                                    SegmentDescriptionByAmira,
+                                    (void *) (segmentIndexList+i)
+                                    )==NULL)
+      return(1);
+
+    /*       sprintf(BndCondName, "BndCond of Segment %d", i); */
+    /*       if(CreateBoundaryCondition(BndCondName, i,DummyBndCondition, NULL) == NULL) */
+    /*       return(1); */
+
+  }
+
+  printf("%d segments created!\n", noOfSegments);
+
+  /* That's it!*/
+  return(0);
+
+
+}
+
+
+int ListEnvCommand (int argc, char **argv);
+
+/** \todo Clear grid before reading! */
+//template<>
+void Dune::AmiraMeshReader<Dune::UGGrid<3,3>, double>::read(Dune::UGGrid<3,3>& grid,
+                                                            const std::string& filename)  try
+{
+  printf("This is the AmiraMesh reader for UGGrid<3,3>!\n");
+
+  //loaddomain $file @PARA_FILE $name @DOMAIN
+  CreateDomain(grid, "olisDomain", "cube.par.am");
+
+  //configure @PROBLEM $d @DOMAIN;
+  char* configureArgs[2] = {"configure DuneDummyProblem", "d olisDomain"};
+  UG3d::ConfigureCommand(2, configureArgs);
+
+
+  printf("new DuneMG $b DuneDummyProblem $f DuneFormat $h 1G\n");
+  //new @PROBLEM $b @PROBLEM $f @FORMAT $h @HEAP;
+  char* newArgs[4] = {"new DuneMG", "b DuneDummyProblem", "f DuneFormat", "h 1G"};
+  if (UG3d::NewCommand(4, newArgs))
+    assert(false);
+
+
+  // ////////////////////////////////////////////
+  // loadmesh $file @GRID_FILE $name @DOMAIN;
+
+  const int DIM = 3;
+
+  int i;
+  double nodePos[DIM];
+#if 0
+  char DomainName[NAMESIZE];
+  char filename[NAMESIZE];
+  INT maxBndNodeID, noOfInnerNodes, noOfBndNodes,
+      noOfElem, noOfCreatedElem, nMapNodes, created;
+#endif
+  UG3d::NODE* theNode;
+
+  printf("Loading Amira mesh %s\n", filename.c_str());
+
+  // /////////////////////////////////////////////////////
+  // Load the AmiraMesh file
+  AmiraMesh* am = AmiraMesh::read(filename.c_str());
+
+  if(!am)
+    throw("Could not open AmiraMesh file");
+
+  float* am_node_coordinates_float = NULL;
+  double* am_node_coordinates_double = NULL;
+
+  // get the different data fields
+  AmiraMesh::Data* am_coordinateData =  am->findData("Nodes", HxFLOAT, 3, "Coordinates");
+  if (am_coordinateData)
+    am_node_coordinates_float = (float*) am_coordinateData->dataPtr();
+  else {
+    am_coordinateData =  am->findData("Nodes", HxDOUBLE, 3, "Coordinates");
+    if (am_coordinateData)
+      am_node_coordinates_double = (double*) am_coordinateData->dataPtr();
+    else
+      throw("No vertex coordinates found in the file!");
+
+  }
+
+
+  AmiraMesh::Data* tetrahedronData = am->findData("Tetrahedra", HxINT32, 4, "Nodes");
+  int*  elemData         = (int*)tetrahedronData->dataPtr();
+
+  /*
+     All Boundary nodes are  assumed to be inserted already.
+     We just have to insert the inner nodes and the elements
+   */
+  UG3d::multigrid* theMG = UG3d::GetMultigrid("DuneMG");
+  assert(theMG);
+
+  int maxBndNodeID = -1;
+  for (theNode=theMG->grids[0]->firstNode[0]; theNode!=NULL; theNode=theNode->succ)
+  {
+    // The following two lines ought to be in here, but the
+    // OBJT macros is somewhat complicated, so I leave it out
+    // for the time being.
+    //       if(OBJT(theNode->myvertex) == UG3d::IVOBJ)
+    //           UserWriteF("Warning: Node %d is inner node\n", ID(theNode));
+    maxBndNodeID = MAX(theNode->id, maxBndNodeID);
+  }
+  UG3d::UserWriteF("Already %d nodes existing\n", maxBndNodeID+1);
+
+
+  //noOfBndNodes = maxBndNodeID;
+
+
+  int noOfNodes = am->nElements("Nodes");
+
+  //  noOfInnerNodes = noOfNodes - noOfBndNodes;
+  printf("AmiraMesh has %d total nodes\n", noOfNodes);
+
+
+
+
+
+  for(i = 0; i < noOfNodes; i++) {
+
+    assert(am_node_coordinates_float || am_node_coordinates_double);
+    if (am_node_coordinates_float) {
+      nodePos[0] = am_node_coordinates_float[3*i];
+      nodePos[1] = am_node_coordinates_float[3*i+1];
+      nodePos[2] = am_node_coordinates_float[3*i+2];
+    } else {
+      nodePos[0] = am_node_coordinates_double[3*i];
+      nodePos[1] = am_node_coordinates_double[3*i+1];
+      nodePos[2] = am_node_coordinates_double[3*i+2];
+    }
+
+    /// \todo Warum ist nicht UG3d::InsertInnerNode Pflicht???
+    if (InsertInnerNode(theMG->grids[0], nodePos) == NULL)
+      throw("inserting an inner node failed");
+
+  }
+
+
+
+  /* all inner nodes are inserted , now we insert the elements */
+  int noOfElem = am->nElements("Tetrahedra");;
+
+  int noOfCreatedElem = 0;
+  for(i=0; i < noOfElem; i++)
+  {
+    int cornerIDs[4], j;
+
+    //AmiraGetCornerIDsOfElem(i, cornerIDs);
+
+    /* only tetrahedrons */
+    /* printf("MeshAccess: elem id : %d, node ids : %d %d %d %d\n", i,
+       elemData[4*i+0], elemData[4*i+1], elemData[4*i+2], elemData[4*i+3]);*/
+    cornerIDs[0] = elemData[4*i]-1;
+    cornerIDs[1] = elemData[4*i+1]-1;
+    cornerIDs[2] = elemData[4*i+2]-1;
+    cornerIDs[3] = elemData[4*i+3]-1;
+
+
+    /*       printf("elem id : %d, node ids : %d %d %d %d\n", i, cornerIDs[0], cornerIDs[1], cornerIDs[2], cornerIDs[3]); */
+
+    if (InsertElementFromIDs(theMG->grids[0], 4,cornerIDs, NULL) == NULL)
+      throw("inserting an element failed");
+
+    noOfCreatedElem++;
+    //       if(statusBar)
+    //         if(i % 150 == 0)
+    //           printf("Already %d of %d elements created: %2.1f %%\r",i,noOfElem,
+    //                      100. * (DOUBLE) i / ((DOUBLE)(noOfElem)));
+
+  }
+
+  //   if(statusBar) {
+  //       UserWriteF("Finished: %d of %d elements created: %2.1f %%\n",i,noOfElem,
+  //                  100. * (DOUBLE) i / ((DOUBLE)(noOfElem)));
+  //       printf("\n");
+  //   }
+
+
+  if(noOfElem != noOfCreatedElem)
+    throw("inserting an element failed");
+
+  UG3d::UserWriteF("amiraloadmesh: %d elements created\n", noOfCreatedElem);
+
+  // set the subdomainIDs
+  AmiraMesh::Data* am_material_ids = am->findData("Tetrahedra", HxBYTE, 1, "Materials");
+  if (!am_material_ids)
+    throw("Field 'Materials' not found.");
+
+  char* material_ids         = (char*)am_material_ids->dataPtr();
+
+  i = 0;
+  UG3d::ELEMENT* theElement;
+  for (theElement=theMG->grids[0]->elements[0]; theElement!=NULL; theElement=theElement->ge.succ)
+  {
+
+    /* get subdomain of element */
+    int id = material_ids[i];
+
+    /** \todo These macros need to go, of course */
+#define ControlWord(p,ce) (((unsigned int *)(p))[UG3d::control_entries[ce].offset_in_object])
+
+#define         CW_WRITE(p, ce, n)   ControlWord(p,ce) = (ControlWord(p,ce)&UG3d::control_entries[ce].xor_mask)|(((n)<<UG3d::control_entries[ce].offset_in_word)&UG3d::control_entries[ce].mask)
+
+#define SETSUBDOMAIN(p,n) CW_WRITE(p,UG3d::SUBDOMAIN_CE,n)
+
+    SETSUBDOMAIN(theElement, id+1);
+
+    //printf("element: %d,  id: %d\n", i, id+1);
+
+    i++;
+  }
+
+  delete am;
+
+
+
+  //   if (MG_COARSE_FIXED(theMG))
+  //     return (GM_OK);
+  UG3d::SetEdgeAndNodeSubdomainFromElements(theMG->grids[0]);
+
+
+  /** \todo Do we really need to call CreateAlgebra for Dune?*/
+  //   if (UG3d::CreateAlgebra(theMG) != GM_OK)
+  //       REP_ERR_RETURN (GM_ERROR);
+
+  /** \todo Check whether this release is necessary */
+  /* here all temp memory since CreateMultiGrid is released */
+  //UG3d::ReleaseTmpMem(MGHEAP(theMG),MG_MARK_KEY(theMG));
+#define FROM_TOP   1
+#define ReleaseTmpMem(p,k) Release(p,FROM_TOP,k)
+  ReleaseTmpMem(theMG->theHeap, theMG->MarkKey);
+  theMG->MarkKey = 0;
+
+  UG3d::UserWriteF("Coarse grid fixed. Do not call fixcoarsegrid.\n");
+
+
+  return;
+
+
+}
+catch (const char* msg) {
+
+  printf("%s\n", msg);
+  return;
+}

lib/io/amuggridwriter.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+// ///////////////////////////////////////////////
+// Specialization of the AmiraMesh writer for SimpleGrid<3,3>
+// ///////////////////////////////////////////////
+
+#include "../../grid/uggrid.hh"
+
+#include <amiramesh/AmiraMesh.h>
+
+namespace Dune {
+
+  template<>
+  class AmiraMeshWriter<UGGrid<3,3>, double> {
+
+  public:
+
+    static void write(const UGGrid<3,3>& grid,
+                      const Array<double>& sol,
+                      const std::string& filename);
+
+
+    AmiraMeshWriter() {}
+
+  };
+
+}
+
+//template<>
+void Dune::AmiraMeshWriter<Dune::UGGrid<3,3>, double>::write(const Dune::UGGrid<3,3>& grid,
+                                                             const Array<double>& sol,
+                                                             const std::string& filename)
+{
+  printf("This is the AmiraMesh writer for UGGrid<3,3>!\n");
+
+  const int DIM = 3;
+
+  int noOfNodes = 0;  //grid.size(0, 3);
+
+  printf("hem\n");
+  int fl, tl, k, i, noOfElem, noOfBndTri, MarkKey, ncomp, maxSubDom;
+  //UG3d::VECDATA_DESC *sol = NULL;
+  UG3d::ELEMENT *t, **elemList;
+  UG3d::VECTOR *vec;
+  UG3d::NODE* theNode;
+  std::string solFilename;
+
+
+  fl = 0;
+  //     tl = TOPLEVEL(theMG);
+
+  // Construct the name for the geometry file
+  std::string geoFilename(filename);
+  geoFilename += ".am";
+
+#if 0
+  /* reset vcflag */
+  for (k=fl; k<=tl; k++)
+    for (vec=PFIRSTVECTOR(GRID_ON_LEVEL(theMG,k)); vec!= NULL; vec=SUCCVC(vec))
+      SETVCFLAG(vec, 0);
+
+
+  /*count nodes and enumerate nodes */
+  noOfNodes = 0;
+  for (k=fl; k<=tl; k++)
+    for (t=FIRSTELEMENT(GRID_ON_LEVEL(theMG,k)); t!=NULL; t=SUCCE(t))
+      if(AmiraElement(t, tl, k))
+        if(EnumerateAmiraNodes(t, sol, &noOfNodes))
+          return(PARAMERRORCODE);
+#endif
+
+#if 0
+  /* build up list of all surface elements and count nodes */
+  noOfElem = 0;
+  maxSubDom = 0;
+  for (k=fl; k<=tl; k++)
+    for (t=FIRSTELEMENT(GRID_ON_LEVEL(theMG,k)); t!=NULL; t=SUCCE(t))
+      if(AmiraElement(t, tl, k))
+      {
+        noOfElem++;
+        maxSubDom = MAX(maxSubDom, SUBDOMAIN(t));
+      }
+
+  MarkTmpMem(MGHEAP(theMG), &MarkKey);
+  elemList = (ELEMENT **) GetTmpMem(MGHEAP(theMG),  (noOfElem+1)*sizeof(ELEMENT *), MarkKey);
+  if(elemList == NULL)
+    return(1);
+
+
+  i = noOfNodes;
+  noOfNodes = 0;
+  for (k=fl; k<=tl; k++)
+    for (vec=PFIRSTVECTOR(GRID_ON_LEVEL(theMG,k)); vec!= NULL; vec=SUCCVC(vec))
+    {
+
+      if(!VCFLAG(vec))
+        continue;
+      VINDEX(vec) = ++noOfNodes;
+      if(VOTYPE(vec) != NODEVEC)
+      {
+        PrintErrorMessage('E',"WriteAmiraGeometry","vector not of type NODEVEC");
+        return (PARAMERRORCODE);
+      }
+    }
+
+  if(FULLREFINELEVEL(theMG) < TOPLEVEL(theMG))
+    for (k = tl; k >= MAX(fl, BOTTOMLEVEL(theMG)+1); k--)
+      for (theNode=FIRSTNODE(GRID_ON_LEVEL(theMG,k)); theNode!=NULL; theNode=SUCCN(theNode))
+      {
+        NODE* fatherNode;
+        VECTOR *fatherVec;
+
+        vec = NVECTOR(theNode);
+
+        if(vec == NULL)
+          continue;
+
+        if(!CORNERTYPE(theNode))
+          continue;
+
+        fatherNode =  (NODE*) NFATHER(theNode);
+
+        if(fatherNode == NULL)
+          continue;
+
+        fatherVec = (VECTOR*) NVECTOR(fatherNode);
+
+        if(fatherVec == NULL)
+        {
+          PrintErrorMessage('E',"WriteAmiraGeometry","father node without vector");
+          return (PARAMERRORCODE);
+        }
+        if(VCFLAG(fatherVec))
+        {
+          PrintErrorMessage('E',"WriteAmiraGeometry","father vec with VCFLAG");
+          return (PARAMERRORCODE);
+        }
+        VINDEX(fatherVec) = VINDEX(vec);
+      }
+
+  if(i != noOfNodes)
+    PrintErrorMessage('E',"WriteAmira","wrong number of nodes");
+
+  assert(i == noOfNodes);
+  /* write elemts to list */
+
+  noOfElem   = 0;
+  noOfBndTri = 0;
+  for (k=fl; k<=tl; k++)
+    for (t=FIRSTELEMENT(GRID_ON_LEVEL(theMG,k)); t!=NULL; t=SUCCE(t))
+      if (AmiraElement(t, tl, k))
+      {
+        elemList[noOfElem++] = t;
+#ifndef ModelP
+        if (OBJT(t) != BEOBJ)
+          continue;
+#endif
+        for (i=0; i<SIDES_OF_ELEM(t); i++)
+          if (ElemSideOnPlotBnd(t, i, NULL) > 0)
+            noOfBndTri++;
+      }
+#endif
+
+
+  // create amiramesh object
+  AmiraMesh am_geometry;
+
+  // write grid vertex coordinates
+  AmiraMesh::Location* geo_nodes = new AmiraMesh::Location("Nodes", noOfNodes);
+  am_geometry.insert(geo_nodes);
+
+  AmiraMesh::Data* geo_node_data = new AmiraMesh::Data("Coordinates", geo_nodes,
+                                                       McPrimType::mc_float, DIM);
+  am_geometry.insert(geo_node_data);
+
+
+  //    i=0;
+  //    for (k=fl; k<=tl; k++)
+  //        for (vec=FIRSTVECTOR(GRID_ON_LEVEL(theMG,k)); vec!= NULL; vec=SUCCVC(vec))
+  //            if ((VOTYPE(vec) == NODEVEC) && ((k==tl) || (VNCLASS(vec)<1) ))
+  //                {
+  //                    /* write coordinates to geo file */
+  //                    DOUBLE pos[DIM];
+
+  //                    if(VectorPosition(vec, pos))
+  //                        return(PARAMERRORCODE);
+
+  //                    ((float*)geo_node_data->dataPtr())[i++] = pos[0];
+  //                    ((float*)geo_node_data->dataPtr())[i++] = pos[1];
+  //                    ((float*)geo_node_data->dataPtr())[i++] = pos[2];
+
+  //                }
+
+  UGGridLevelIterator<3, 3, 3> vertex = grid.lbegin<3>(0);
+  printf("level: %d\n", vertex.level());
+
+
+
+  UGGridLevelIterator<0, 3, 3> element = grid.lbegin<0>(0);
+
+  for (i=0; element!=grid.lend<0>(0); element++, i++)
+    printf("Element! %d\n", i);
+
+#if 0
+  // handle materials
+  HxParamBundle* materials = new HxParamBundle("Materials");
+
+  for (k=0; k<=maxSubDom; k++) {
+
+    char buffer[100];
+    sprintf(buffer, "Material%d", k);
+    HxParamBundle* newMaterial = new HxParamBundle(buffer);
+
+    HxParameter* newId = new HxParameter("Id", k);
+    newMaterial->insert(newId);
+
+    materials->insert(newMaterial);
+
+  }
+
+  am_geometry.parameters.insert(materials);
+
+  ncomp = 0;
+  for(i=0; i<NVECTYPES; i++)
+    ncomp = MAX(ncomp,VD_NCMPS_IN_TYPE(sol, i));
+
+
+  /* write element section to geo - file */
+  AmiraMesh::Location* element_loc = new AmiraMesh::Location("Tetrahedra", noOfElem);
+  am_geometry.insert(element_loc);
+
+  AmiraMesh::Data* element_data = new AmiraMesh::Data("Nodes", element_loc,
+                                                      McPrimType::mc_int32, DIM+1);
+  am_geometry.insert(element_data);
+
+  int (*dPtr)[DIM+1] = (int(*)[DIM+1])element_data->dataPtr();
+
+  for(i=0; i<noOfElem; i++)
+    if(WriteAmiraGeometry(dPtr[i], elemList[i], sol))
+      return(PARAMERRORCODE);
+
+  // write material section to geo-file
+  AmiraMesh::Data* element_materials = new AmiraMesh::Data("Materials", element_loc, McPrimType::mc_uint8, 1);
+  am_geometry.insert(element_materials);
+
+  for(i=0; i<noOfElem; i++)
+    ((unsigned char*)element_materials->dataPtr())[i] = SUBDOMAIN(elemList[i]);
+
+  //    fprintf(geoFile, "\n@%d\n", 4 /* BND_TRIDATA_SECTION */);
+  //    for(i=0;i<noOfElem;i++)
+  //      if(WriteAmiraBndTriData(geoFile, elemList[i], ncomp))
+  //        return(PARAMERRORCODE);
+
+  //    fprintf(geoFile, "\n@%d\n", 5 /* BND_TRI_SECTION */);
+  //    for(i=0;i<noOfElem;i++)
+  //      if(WriteAmiraBndTri(geoFile, elemList[i]))
+  //        return(PARAMERRORCODE);
+
+
+
+
+  //////////////////////////////////////////////////////
+
+  // Now save the solution
+  // Construct the name for the solution file
+
+  //strcpy(solFilename, (const char*) BaseName);
+#ifdef ModelP
+  //strcat(solFilename, procNumber);
+#endif
+  //strcat(solFilename, ".sol");
+
+  //AmiraMesh am_solution;
+
+  //AmiraMesh::Location* nodes = new AmiraMesh::Location("Nodes", noOfNodes);
+  //am_solution.insert(nodes);
+
+  AmiraMesh::Data* nodeData = new AmiraMesh::Data("Data", geo_nodes, McPrimType::mc_float, ncomp);
+  am_geometry.insert(nodeData);
+
+  AmiraMesh::Field* nodeField = new AmiraMesh::Field("f", ncomp, McPrimType::mc_float,
+                                                     AmiraMesh::t_linear, nodeData);
+  am_geometry.insert(nodeField);
+
+  i=0;
+  for (k=fl; k<=tl; k++)
+    for (vec=FIRSTVECTOR(GRID_ON_LEVEL(theMG,k)); vec!= NULL; vec=SUCCVC(vec))
+      if ((VOTYPE(vec) == NODEVEC) && ((k==tl) || (VNCLASS(vec)<1) ))
+      {
+        SHORT vtype   = VTYPE(vec);
+        /** \bug This definition of ncomp shadows another one! */
+        SHORT ncomp   = VD_NCMPS_IN_TYPE(sol, vtype);
+        SHORT *cmpptr = VD_CMPPTR_OF_TYPE(sol, vtype);
+        SHORT j;
+
+        for (j=0; j<ncomp; j++)
+          ((float*)nodeData->dataPtr())[ncomp*i+j] = VVALUE(vec, cmpptr[j]);
+
+        i++;
+      }
+
+  // actually save the solution file
+  //    if (!am_solution.write(solFilename, 1) ) {
+  //        printf("An error has occured writing file %s.\n", solFilename);
+  //        return PARAMERRORCODE;
+  //    }
+
+#endif
+  if(!am_geometry.write(geoFilename.c_str(), 1)) {
+    printf("writing geometry file failed in amira : \n");
+    /** \todo Do a decent error handling */
+    return;
+  }
+
+  printf("Grid written successfully to:\n %s \n", geoFilename.c_str());
+}