From 946b8b3f4b5bb26091a4a557af61a240ab6c1185 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@dune-project.org>
Date: Tue, 2 Mar 2004 11:15:05 +0000
Subject: [PATCH] I/O for the amiramesh format; no functionality so far
[[Imported from SVN: r357]]
---
lib/io/amirameshreader.hh | 50 ++++
lib/io/amirameshwriter.hh | 9 +
lib/io/amuggridreader.cc | 470 ++++++++++++++++++++++++++++++++++++++
lib/io/amuggridwriter.cc | 316 +++++++++++++++++++++++++
4 files changed, 845 insertions(+)
create mode 100644 lib/io/amirameshreader.hh
create mode 100644 lib/io/amuggridreader.cc
create mode 100644 lib/io/amuggridwriter.cc
diff --git a/lib/io/amirameshreader.hh b/lib/io/amirameshreader.hh
new file mode 100644
index 000000000..111b7fe73
--- /dev/null
+++ b/lib/io/amirameshreader.hh
@@ -0,0 +1,50 @@
+// -*- 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
diff --git a/lib/io/amirameshwriter.hh b/lib/io/amirameshwriter.hh
index 8e77764f5..a51c8bd2a 100644
--- a/lib/io/amirameshwriter.hh
+++ b/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
diff --git a/lib/io/amuggridreader.cc b/lib/io/amuggridreader.cc
new file mode 100644
index 000000000..8bff05516
--- /dev/null
+++ b/lib/io/amuggridreader.cc
@@ -0,0 +1,470 @@
+// -*- 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;
+}
diff --git a/lib/io/amuggridwriter.cc b/lib/io/amuggridwriter.cc
new file mode 100644
index 000000000..bd6e4d50d
--- /dev/null
+++ b/lib/io/amuggridwriter.cc
@@ -0,0 +1,316 @@
+// -*- 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());
+}
--
GitLab