src/escriptreader/MeshWithElements.cpp


/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/

//
// MeshWithElements.cpp
//
#include <escriptreader/MeshWithElements.h>
#include <escriptreader/ElementData.h>
#include <netcdf.hh>
#if HAVE_SILO
#include <silo.h>
#endif

using namespace std;

namespace EscriptReader {

//
//
//
MeshWithElements::MeshWithElements() : Mesh(),
    cells(NULL),
    faces(NULL),
    contacts(NULL),
    points(NULL)
{
    name = "Elements";
}

//
//
//
MeshWithElements::MeshWithElements(const MeshWithElements& m) :
    Mesh(m)
{
    nodeTag = m.nodeTag;
    nodeGDOF = m.nodeGDOF;
    nodeGNI = m.nodeGNI;
    nodeGRDFI = m.nodeGRDFI;
    nodeGRNI = m.nodeGRNI;
    cells = new ElementData(*m.cells);
    faces = new ElementData(*m.faces);
    contacts = new ElementData(*m.contacts);
    points = new ElementData(*m.points);
}

//
//
//
MeshWithElements::~MeshWithElements()
{
    delete cells;
    delete faces;
    delete contacts;
    delete points;
}

//
// Returns a vector of strings containing Silo mesh names that have been
// written
//
StringVec MeshWithElements::getMeshNames() const
{
    StringVec res;
    res.push_back(name);
    StringVec tmpVec;
    tmpVec = cells->getMeshNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = faces->getMeshNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = contacts->getMeshNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = points->getMeshNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    return res;
}

Mesh* MeshWithElements::getMeshByName(const string name) const
{
    Mesh* ret = NULL;
    if (name == "Elements")
        ret = cells->fullMesh;
    else if (name == "ReducedElements")
        ret = cells->reducedMesh;
    else if (name == "FaceElements")
        ret = faces->fullMesh;
    else if (name == "ReducedFaceElements")
        ret = faces->reducedMesh;
    else if (name == "ContactElements")
        ret = contacts->fullMesh;
    else if (name == "ReducedContactElements")
        ret = contacts->reducedMesh;
    else if (name == "Points")
        ret = points->fullMesh;

    return ret;
}

//
// Returns a vector of strings containing Silo variable names that have
// been written
//
StringVec MeshWithElements::getVarNames() const
{
    StringVec res;
    res.push_back("Nodes_Id");
    res.push_back("Nodes_Tag");
    res.push_back("Nodes_gDOF");
    res.push_back("Nodes_gNI");
    res.push_back("Nodes_grDfI");
    res.push_back("Nodes_grNI");

    StringVec tmpVec;
    tmpVec = cells->getVarNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = faces->getVarNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = contacts->getVarNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());
    tmpVec = points->getVarNames();
    res.insert(res.end(), tmpVec.begin(), tmpVec.end());

    return res;
}

const IntVec& MeshWithElements::getVarDataByName(const string name) const
{
    if (name == "Nodes_Id")
        return nodeID;
    else if (name == "Nodes_Tag")
        return nodeTag;
    else if (name == "Nodes_gDOF")
        return nodeGDOF;
    else if (name == "Nodes_gNI")
        return nodeGNI;
    else if (name == "Nodes_grDfI")
        return nodeGRDFI;
    else if (name == "Nodes_grNI")
        return nodeGRNI;
    else if (name.compare(0, 9, "Elements_") == 0)
        return cells->getVarDataByName(name);
    else if (name.compare(0, 13, "FaceElements_") == 0)
        return faces->getVarDataByName(name);
    else if (name.compare(0, 16, "ContactElements_") == 0)
        return contacts->getVarDataByName(name);
    else if (name.compare(0, 7, "Points_") == 0)
        return points->getVarDataByName(name);
    else
        throw "Invalid variable name";
}

ElementData* MeshWithElements::getElementsByName(const string name) const
{
    ElementData* ret = NULL;
    if (name == "Elements" || name == "ReducedElements")
        ret = cells;
    else if (name == "FaceElements" || name == "ReducedFaceElements")
        ret = faces;
    else if (name == "ContactElements" || name == "ReducedContactElements")
        ret = contacts;
    else if (name == "Points")
        ret = points;

    return ret;
}

//
// Reads mesh and element data from NetCDF file with given name
//
bool MeshWithElements::readFromNc(const string& filename)
{
    if (!Mesh::readFromNc(filename))
        return false;

    NcError ncerr(NcError::silent_nonfatal);
    NcFile* input;
    NcVar* var;
 
    input = new NcFile(filename.c_str());
    if (!input->is_valid()) {
        cerr << "Could not open input file " << filename << "." << endl;
        delete input;
        return false;
    }

    nodeTag.clear();
    nodeTag.insert(nodeTag.end(), numNodes, 0);
    var = input->get_var("Nodes_Tag");
    var->get(&nodeTag[0], numNodes);

    nodeGDOF.clear();
    nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
    var = input->get_var("Nodes_gDOF");
    var->get(&nodeGDOF[0], numNodes);

    nodeGNI.clear();
    nodeGNI.insert(nodeGNI.end(), numNodes, 0);
    var = input->get_var("Nodes_gNI");
    var->get(&nodeGNI[0], numNodes);

    nodeGRDFI.clear();
    nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
    var = input->get_var("Nodes_grDfI");
    var->get(&nodeGRDFI[0], numNodes);

    nodeGRNI.clear();
    nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
    var = input->get_var("Nodes_grNI");
    var->get(&nodeGRNI[0], numNodes);

    // Read all element types
    cells = new ElementData("Elements", this);
    cells->readFromNc(input);
    faces = new ElementData("FaceElements", this);
    faces->readFromNc(input);
    contacts = new ElementData("ContactElements", this);
    contacts->readFromNc(input);
    points = new ElementData("Points", this);
    points->readFromNc(input);

    delete input;
    return true;
}

//
//
//
void MeshWithElements::handleGhostZones(int ownIndex)
{
    cells->handleGhostZones(ownIndex);
    faces->handleGhostZones(ownIndex);
    contacts->handleGhostZones(ownIndex);
    points->handleGhostZones(ownIndex);
#ifdef _DEBUG
    cout << "block " << ownIndex << " has " << cells->getGhostCount()
         << " (" << cells->getReducedGhostCount() << ") ghost zones," << endl;
    cout << faces->getGhostCount() << " (" << faces->getReducedGhostCount()
         << ") ghost faces," << endl;
    cout << contacts->getGhostCount() << " ("
         << contacts->getReducedGhostCount() << ") ghost contacts," << endl;
    cout << points->getGhostCount() << " (" << points->getReducedGhostCount()
         << ") ghost points." << endl;
#endif
}

//
//
//
void MeshWithElements::removeGhostZones()
{
    cells->removeGhostZones();
    faces->removeGhostZones();
    contacts->removeGhostZones();
    points->removeGhostZones();
#ifdef _DEBUG
    cout << "After removing ghost zones there are" << endl;
    cout << "    " << numNodes << " Nodes, ";
    cout << cells->getCount() << " Elements, ";
    cout << faces->getCount() << " Face elements, ";
    cout << contacts->getCount() << " Contact elements, ";
    cout << points->getCount() << " Points left." << endl;
#endif
}

//
//
//
bool MeshWithElements::writeToSilo(DBfile* dbfile, const string& pathInSilo)
{
#if HAVE_SILO
    int ret;

    // Write meshes and zone-centered variables
    if (!cells->writeToSilo(dbfile, pathInSilo) ||
            !faces->writeToSilo(dbfile, pathInSilo) ||
            !contacts->writeToSilo(dbfile, pathInSilo) ||
            !points->writeToSilo(dbfile, pathInSilo))
        return false;

    if (!Mesh::writeToSilo(dbfile, pathInSilo))
        return false;
    
    if (pathInSilo != "") {
        ret = DBSetDir(dbfile, pathInSilo.c_str());
        if (ret != 0)
            return false;
    }

    string siloMeshName = getFullSiloName();

    // Write node-centered variables
    ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
                (float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);
    if (ret == 0)
        ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
                (float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);
    if (ret == 0)
        ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
                (float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);
    if (ret == 0)
        ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
                (float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);
    if (ret == 0)
        ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
                (float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
                DB_NODECENT, NULL);

    DBSetDir(dbfile, "/");
    return (ret == 0);

#else // !HAVE_SILO
    return false;
#endif
}

} // namespace EscriptReader

1	caltinay	2183
2			/*******************************************************
3			*
4			* Copyright (c) 2003-2008 by University of Queensland
5			* Earth Systems Science Computational Center (ESSCC)
6			* http://www.uq.edu.au/esscc
7			*
8			* Primary Business: Queensland, Australia
9			* Licensed under the Open Software License version 3.0
10			* http://www.opensource.org/licenses/osl-3.0.php
11			*
12			*******************************************************/
13
14			//
15			// MeshWithElements.cpp
16			//
17			#include <escriptreader/MeshWithElements.h>
18			#include <escriptreader/ElementData.h>
19			#include <netcdf.hh>
20			#if HAVE_SILO
21			#include <silo.h>
22			#endif
23
24			using namespace std;
25
26	caltinay	2187	namespace EscriptReader {
27
28	caltinay	2183	//
29			//
30			//
31	caltinay	2194	MeshWithElements::MeshWithElements() : Mesh(),
32			cells(NULL),
33			faces(NULL),
34			contacts(NULL),
35			points(NULL)
36	caltinay	2183	{
37			name = "Elements";
38			}
39
40			//
41			//
42			//
43			MeshWithElements::MeshWithElements(const MeshWithElements& m) :
44			Mesh(m)
45			{
46			nodeTag = m.nodeTag;
47			nodeGDOF = m.nodeGDOF;
48			nodeGNI = m.nodeGNI;
49			nodeGRDFI = m.nodeGRDFI;
50			nodeGRNI = m.nodeGRNI;
51			cells = new ElementData(*m.cells);
52			faces = new ElementData(*m.faces);
53			contacts = new ElementData(*m.contacts);
54			points = new ElementData(*m.points);
55			}
56
57			//
58			//
59			//
60			MeshWithElements::~MeshWithElements()
61			{
62			delete cells;
63			delete faces;
64			delete contacts;
65			delete points;
66			}
67
68			//
69			// Returns a vector of strings containing Silo mesh names that have been
70			// written
71			//
72			StringVec MeshWithElements::getMeshNames() const
73			{
74			StringVec res;
75			res.push_back(name);
76			StringVec tmpVec;
77			tmpVec = cells->getMeshNames();
78			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
79			tmpVec = faces->getMeshNames();
80			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
81			tmpVec = contacts->getMeshNames();
82			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
83			tmpVec = points->getMeshNames();
84			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
85			return res;
86			}
87
88			Mesh* MeshWithElements::getMeshByName(const string name) const
89			{
90			Mesh* ret = NULL;
91			if (name == "Elements")
92			ret = cells->fullMesh;
93			else if (name == "ReducedElements")
94			ret = cells->reducedMesh;
95			else if (name == "FaceElements")
96			ret = faces->fullMesh;
97			else if (name == "ReducedFaceElements")
98			ret = faces->reducedMesh;
99			else if (name == "ContactElements")
100			ret = contacts->fullMesh;
101			else if (name == "ReducedContactElements")
102			ret = contacts->reducedMesh;
103			else if (name == "Points")
104			ret = points->fullMesh;
105
106			return ret;
107			}
108
109			//
110			// Returns a vector of strings containing Silo variable names that have
111			// been written
112			//
113			StringVec MeshWithElements::getVarNames() const
114			{
115			StringVec res;
116			res.push_back("Nodes_Id");
117			res.push_back("Nodes_Tag");
118			res.push_back("Nodes_gDOF");
119			res.push_back("Nodes_gNI");
120			res.push_back("Nodes_grDfI");
121			res.push_back("Nodes_grNI");
122
123			StringVec tmpVec;
124			tmpVec = cells->getVarNames();
125			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
126			tmpVec = faces->getVarNames();
127			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
128			tmpVec = contacts->getVarNames();
129			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
130			tmpVec = points->getVarNames();
131			res.insert(res.end(), tmpVec.begin(), tmpVec.end());
132
133			return res;
134			}
135
136			const IntVec& MeshWithElements::getVarDataByName(const string name) const
137			{
138			if (name == "Nodes_Id")
139			return nodeID;
140			else if (name == "Nodes_Tag")
141			return nodeTag;
142			else if (name == "Nodes_gDOF")
143			return nodeGDOF;
144			else if (name == "Nodes_gNI")
145			return nodeGNI;
146			else if (name == "Nodes_grDfI")
147			return nodeGRDFI;
148			else if (name == "Nodes_grNI")
149			return nodeGRNI;
150			else if (name.compare(0, 9, "Elements_") == 0)
151			return cells->getVarDataByName(name);
152			else if (name.compare(0, 13, "FaceElements_") == 0)
153			return faces->getVarDataByName(name);
154			else if (name.compare(0, 16, "ContactElements_") == 0)
155			return contacts->getVarDataByName(name);
156			else if (name.compare(0, 7, "Points_") == 0)
157			return points->getVarDataByName(name);
158			else
159	caltinay	2201	throw "Invalid variable name";
160	caltinay	2183	}
161
162			ElementData* MeshWithElements::getElementsByName(const string name) const
163			{
164			ElementData* ret = NULL;
165			if (name == "Elements" \|\| name == "ReducedElements")
166			ret = cells;
167			else if (name == "FaceElements" \|\| name == "ReducedFaceElements")
168			ret = faces;
169			else if (name == "ContactElements" \|\| name == "ReducedContactElements")
170			ret = contacts;
171			else if (name == "Points")
172			ret = points;
173
174			return ret;
175			}
176
177			//
178			// Reads mesh and element data from NetCDF file with given name
179			//
180			bool MeshWithElements::readFromNc(const string& filename)
181			{
182			if (!Mesh::readFromNc(filename))
183			return false;
184
185			NcError ncerr(NcError::silent_nonfatal);
186			NcFile* input;
187			NcVar* var;
188
189			input = new NcFile(filename.c_str());
190			if (!input->is_valid()) {
191			cerr << "Could not open input file " << filename << "." << endl;
192			delete input;
193			return false;
194			}
195
196			nodeTag.clear();
197			nodeTag.insert(nodeTag.end(), numNodes, 0);
198			var = input->get_var("Nodes_Tag");
199			var->get(&nodeTag[0], numNodes);
200
201			nodeGDOF.clear();
202			nodeGDOF.insert(nodeGDOF.end(), numNodes, 0);
203			var = input->get_var("Nodes_gDOF");
204			var->get(&nodeGDOF[0], numNodes);
205
206			nodeGNI.clear();
207			nodeGNI.insert(nodeGNI.end(), numNodes, 0);
208			var = input->get_var("Nodes_gNI");
209			var->get(&nodeGNI[0], numNodes);
210
211			nodeGRDFI.clear();
212			nodeGRDFI.insert(nodeGRDFI.end(), numNodes, 0);
213			var = input->get_var("Nodes_grDfI");
214			var->get(&nodeGRDFI[0], numNodes);
215
216			nodeGRNI.clear();
217			nodeGRNI.insert(nodeGRNI.end(), numNodes, 0);
218			var = input->get_var("Nodes_grNI");
219			var->get(&nodeGRNI[0], numNodes);
220
221			// Read all element types
222			cells = new ElementData("Elements", this);
223			cells->readFromNc(input);
224			faces = new ElementData("FaceElements", this);
225			faces->readFromNc(input);
226			contacts = new ElementData("ContactElements", this);
227			contacts->readFromNc(input);
228			points = new ElementData("Points", this);
229			points->readFromNc(input);
230
231			delete input;
232			return true;
233			}
234
235			//
236			//
237			//
238			void MeshWithElements::handleGhostZones(int ownIndex)
239			{
240			cells->handleGhostZones(ownIndex);
241			faces->handleGhostZones(ownIndex);
242			contacts->handleGhostZones(ownIndex);
243			points->handleGhostZones(ownIndex);
244			#ifdef _DEBUG
245			cout << "block " << ownIndex << " has " << cells->getGhostCount()
246			<< " (" << cells->getReducedGhostCount() << ") ghost zones," << endl;
247			cout << faces->getGhostCount() << " (" << faces->getReducedGhostCount()
248			<< ") ghost faces," << endl;
249			cout << contacts->getGhostCount() << " ("
250			<< contacts->getReducedGhostCount() << ") ghost contacts," << endl;
251			cout << points->getGhostCount() << " (" << points->getReducedGhostCount()
252			<< ") ghost points." << endl;
253			#endif
254			}
255
256			//
257			//
258			//
259			void MeshWithElements::removeGhostZones()
260			{
261			cells->removeGhostZones();
262			faces->removeGhostZones();
263			contacts->removeGhostZones();
264			points->removeGhostZones();
265			#ifdef _DEBUG
266			cout << "After removing ghost zones there are" << endl;
267			cout << " " << numNodes << " Nodes, ";
268			cout << cells->getCount() << " Elements, ";
269			cout << faces->getCount() << " Face elements, ";
270			cout << contacts->getCount() << " Contact elements, ";
271			cout << points->getCount() << " Points left." << endl;
272			#endif
273			}
274
275			//
276			//
277			//
278			bool MeshWithElements::writeToSilo(DBfile* dbfile, const string& pathInSilo)
279			{
280			#if HAVE_SILO
281			int ret;
282
283			// Write meshes and zone-centered variables
284			if (!cells->writeToSilo(dbfile, pathInSilo) \|\|
285			!faces->writeToSilo(dbfile, pathInSilo) \|\|
286			!contacts->writeToSilo(dbfile, pathInSilo) \|\|
287			!points->writeToSilo(dbfile, pathInSilo))
288			return false;
289
290			if (!Mesh::writeToSilo(dbfile, pathInSilo))
291			return false;
292
293			if (pathInSilo != "") {
294			ret = DBSetDir(dbfile, pathInSilo.c_str());
295			if (ret != 0)
296			return false;
297			}
298
299			string siloMeshName = getFullSiloName();
300
301			// Write node-centered variables
302			ret = DBPutUcdvar1(dbfile, "Nodes_Tag", siloMeshName.c_str(),
303			(float*)&nodeTag[0], numNodes, NULL, 0, DB_INT,
304			DB_NODECENT, NULL);
305			if (ret == 0)
306			ret = DBPutUcdvar1(dbfile, "Nodes_gDOF", siloMeshName.c_str(),
307			(float*)&nodeGDOF[0], numNodes, NULL, 0, DB_INT,
308			DB_NODECENT, NULL);
309			if (ret == 0)
310			ret = DBPutUcdvar1(dbfile, "Nodes_gNI", siloMeshName.c_str(),
311			(float*)&nodeGNI[0], numNodes, NULL, 0, DB_INT,
312			DB_NODECENT, NULL);
313			if (ret == 0)
314			ret = DBPutUcdvar1(dbfile, "Nodes_grDfI", siloMeshName.c_str(),
315			(float*)&nodeGRDFI[0], numNodes, NULL, 0, DB_INT,
316			DB_NODECENT, NULL);
317			if (ret == 0)
318			ret = DBPutUcdvar1(dbfile, "Nodes_grNI", siloMeshName.c_str(),
319			(float*)&nodeGRNI[0], numNodes, NULL, 0, DB_INT,
320			DB_NODECENT, NULL);
321
322			DBSetDir(dbfile, "/");
323			return (ret == 0);
324
325			#else // !HAVE_SILO
326			return false;
327			#endif
328			}
329
330	caltinay	2187	} // namespace EscriptReader
331