/[escript]/branches/trilinos_from_5897/dudley/src/Mesh_write.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/Mesh_write.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 6078 by caltinay, Wed Mar 2 04:13:26 2016 UTC revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 16  Line 16 
16    
17  #include "Mesh.h"  #include "Mesh.h"
18    
19    #include <iomanip>
20    
21    using std::cout;
22    using std::endl;
23    using std::ios;
24    using std::setw;
25    using std::string;
26    
27  namespace dudley {  namespace dudley {
28    
29  /// writes the mesh to the external file fname using the Dudley file format  // private
30  void Dudley_Mesh_write(Dudley_Mesh* in, char* fname)  void Mesh::writeElementInfo(std::ostream& stream, const ElementFile* e,
31                                const string& defaultType) const
32  {  {
33      char error_msg[1024];      if (e != NULL) {
34      FILE *f;          stream << e->ename << " " << e->numElements << endl;
35      int NN, i, j, numDim;          const int NN = e->numNodes;
36      Dudley_TagMap *tag_map = in->TagMap;          for (index_t i=0; i < e->numElements; i++) {
37                stream << e->Id[i] << " " << e->Tag[i];
38      if (in->MPIInfo->size > 1)              for (int j=0; j<NN; j++)
39          throw DudleyException("Mesh_write: only single processor runs are supported.");                  stream << " " << Nodes->Id[e->Nodes[INDEX2(j,i,NN)]];
40                stream << endl;
     /* open file */  
     f = fopen(fname, "w");  
     if (f == NULL)  
     {  
         sprintf(error_msg, "Mesh_write: Opening file %s for writing failed.", fname);  
         throw DudleyException(error_msg);  
     }  
   
     /* write header */  
   
     fprintf(f, "%s\n", in->Name);  
   
     /*  write nodes: */  
     if (in->Nodes != NULL) {  
         numDim = Dudley_Mesh_getDim(in);  
         fprintf(f, "%1dD-Nodes %d\n", numDim, in->Nodes->numNodes);  
         for (i = 0; i < in->Nodes->numNodes; i++)  
         {  
             fprintf(f, "%d %d %d", in->Nodes->Id[i], in->Nodes->globalDegreesOfFreedom[i], in->Nodes->Tag[i]);  
             for (j = 0; j < numDim; j++)  
                 fprintf(f, " %20.15e", in->Nodes->Coordinates[INDEX2(j, i, numDim)]);  
             fprintf(f, "\n");  
41          }          }
42      } else {      } else {
43          fprintf(f, "0D-Nodes 0\n");          stream << defaultType << " 0" << endl;
44      }      }
45    }
46    
47      /*  write elements: */  // private
48      if (in->Elements != NULL) {  void Mesh::printElementInfo(const ElementFile* e, const string& title,
49          fprintf(f, "%s %d\n", in->Elements->ename /*referenceElementSet->referenceElement->Type->Name */ ,                              const string& defaultType, bool full) const
50                  in->Elements->numElements);  {
51          NN = in->Elements->numNodes;      if (e != NULL) {
52          for (i = 0; i < in->Elements->numElements; i++) {          dim_t mine=0, overlap=0;
53              fprintf(f, "%d %d", in->Elements->Id[i], in->Elements->Tag[i]);          for (index_t i=0; i < e->numElements; i++) {
54              for (j = 0; j < NN; j++)              if (e->Owner[i] == MPIInfo->rank)
55                  fprintf(f, " %d", in->Nodes->Id[in->Elements->Nodes[INDEX2(j, i, NN)]]);                  mine++;
56              fprintf(f, "\n");              else
57                    overlap++;
58            }
59            cout << "\t" << title << ": "
60                << e->ename << " " << e->numElements << " (TypeId=" << e->etype
61                << ") owner=" << mine << " overlap=" << overlap << endl;
62            if (full) {
63                const int NN = e->numNodes;
64                cout << "\t     Id   Tag Owner Color:  Nodes" << endl;
65                for (index_t i=0; i < e->numElements; i++) {
66                    cout << "\t" << setw(7) << e->Id[i]
67                         << setw(6) << e->Tag[i]
68                         << setw(6) << e->Owner[i]
69                         << setw(6) << e->Color[i] << ": ";
70                    for (int j=0; j<NN; j++)
71                        cout << setw(6) << Nodes->Id[e->Nodes[INDEX2(j,i,NN)]];
72                    cout << endl;
73                }
74          }          }
75      } else {      } else {
76          fprintf(f, "Tet4 0\n");          cout << "\t" << title << ": " << defaultType << " 0" << endl;
77      }      }
78    }
79    
80      /*  write face elements: */  
81      if (in->FaceElements != NULL) {  void Mesh::write(const std::string& filename) const
82          fprintf(f, "%s %d\n", in->FaceElements->ename /*referenceElementSet->referenceElement->Type->Name */ ,  {
83                  in->FaceElements->numElements);      if (MPIInfo->size > 1)
84          NN = in->FaceElements->numNodes;          throw escript::NotImplementedError("Mesh::write: only single rank "
85          for (i = 0; i < in->FaceElements->numElements; i++) {                                             "runs are supported.");
86              fprintf(f, "%d %d", in->FaceElements->Id[i], in->FaceElements->Tag[i]);  
87              for (j = 0; j < NN; j++)      std::ofstream f(filename.c_str());
88                  fprintf(f, " %d", in->Nodes->Id[in->FaceElements->Nodes[INDEX2(j, i, NN)]]);      if (!f.is_open()) {
89              fprintf(f, "\n");          std::stringstream ss;
90          }          ss << "Mesh::write: Opening file " << filename << " for writing failed";
91      } else {          throw escript::IOError(ss.str());
         fprintf(f, "Tri3 0\n");  
92      }      }
93    
94      /*  write points: */      // write header
95      if (in->Points != NULL) {      f << m_name << endl;
96          fprintf(f, "%s %d\n", in->Points->ename /*referenceElementSet->referenceElement->Type->Name */ ,  
97                  in->Points->numElements);      //  write nodes
98          for (i = 0; i < in->Points->numElements; i++) {      if (Nodes != NULL) {
99              fprintf(f, "%d %d %d\n", in->Points->Id[i], in->Points->Tag[i],          const int numDim = getDim();
100                      in->Nodes->Id[in->Points->Nodes[INDEX2(0, i, 1)]]);          f << numDim << "D-Nodes " << Nodes->getNumNodes() << endl;
101            for (index_t i=0; i<Nodes->getNumNodes(); i++) {
102                f << Nodes->Id[i] << " " << Nodes->globalDegreesOfFreedom[i]
103                  << " " << Nodes->Tag[i];
104                f.setf(ios::scientific, ios::floatfield);
105                f.precision(15);
106                for (int j=0; j<numDim; j++)
107                    f << " " << Nodes->Coordinates[INDEX2(j,i,numDim)];
108                f << endl;
109          }          }
110      } else {      } else {
111          fprintf(f, "Point1 0\n");          f << "0D-Nodes 0" << endl;
112      }      }
113    
114      /*  write tags: */      // write elements
115      if (tag_map) {      writeElementInfo(f, Elements, "Tet4");
116          fprintf(f, "Tags\n");  
117          while (tag_map) {      // write face elements
118              fprintf(f, "%s %d\n", tag_map->name, tag_map->tag_key);      writeElementInfo(f, FaceElements, "Tri3");
119              tag_map = tag_map->next;  
120        // write points
121        writeElementInfo(f, Points, "Point1");
122    
123        // write tags
124        if (tagMap.size() > 0) {
125            f <<  "Tags" << endl;
126            TagMap::const_iterator it;
127            for (it=tagMap.begin(); it!=tagMap.end(); it++) {
128                f << it->first << " " << it->second << endl;
129          }          }
130      }      }
131      fclose(f);      f.close();
132  #ifdef Dudley_TRACE  #ifdef Dudley_TRACE
133      printf("mesh %s has been written to file %s\n", in->Name, fname);      cout << "mesh " << m_name << " has been written to file " << filename << endl;
134  #endif  #endif
135  }  }
136    
137  void Dudley_PrintMesh_Info(Dudley_Mesh * in, bool full)  void Mesh::printInfo(bool full)
138  {  {
139      int NN, i, j, numDim;      cout << "PrintMeshInfo running on CPU " << MPIInfo->rank << " of "
140      Dudley_TagMap *tag_map = in->TagMap;                << MPIInfo->size << endl;
141        cout << "\tMesh name '" << m_name << "'\n";
142      fprintf(stdout, "Dudley_PrintMesh_Info running on CPU %d of %d\n", in->MPIInfo->rank, in->MPIInfo->size);      cout << "\tApproximation order " << approximationOrder << endl;
143      fprintf(stdout, "\tMesh name '%s'\n", in->Name);      cout << "\tIntegration order " << integrationOrder << endl;
144      fprintf(stdout, "\tApproximation order %d\n", in->approximationOrder);      cout << "\tReduced Integration order " << reducedIntegrationOrder << endl;
145      fprintf(stdout, "\tReduced Approximation order %d\n", in->reducedApproximationOrder);  
146      fprintf(stdout, "\tIntegration order %d\n", in->integrationOrder);      // write nodes
147      fprintf(stdout, "\tReduced Integration order %d\n", in->reducedIntegrationOrder);      if (Nodes != NULL) {
148            const int numDim = getDim();
149      /* write nodes: */          cout << "\tNodes: " << numDim << "D-Nodes " << Nodes->getNumNodes() << endl;
     if (in->Nodes != NULL) {  
         numDim = Dudley_Mesh_getDim(in);  
         fprintf(stdout, "\tNodes: %1dD-Nodes %d\n", numDim, in->Nodes->numNodes);  
150          if (full) {          if (full) {
151              fprintf(stdout, "\t     Id   Tag  gDOF   gNI grDfI  grNI:  Coordinates\n");              cout << "\t     Id   Tag  gDOF   gNI grDfI  grNI:  Coordinates\n";
152              for (i = 0; i < in->Nodes->numNodes; i++) {              for (index_t i=0; i < Nodes->getNumNodes(); i++) {
153                  fprintf(stdout, "\t  %5d %5d %5d %5d %5d %5d: ", in->Nodes->Id[i], in->Nodes->Tag[i],                  cout << "\t" << setw(7) << Nodes->Id[i]
154                          in->Nodes->globalDegreesOfFreedom[i], in->Nodes->globalNodesIndex[i],                       << setw(6) << Nodes->Tag[i]
155                          in->Nodes->globalReducedDOFIndex[i], in->Nodes->globalReducedNodesIndex[i]);                       << setw(6) << Nodes->globalDegreesOfFreedom[i]
156                  for (j = 0; j < numDim; j++)                       << setw(6) << Nodes->globalNodesIndex[i]
157                      fprintf(stdout, " %20.15e", in->Nodes->Coordinates[INDEX2(j, i, numDim)]);                       << setw(6) << Nodes->globalDegreesOfFreedom[i]
158                  fprintf(stdout, "\n");                       << setw(6) << Nodes->globalNodesIndex[i] << ": ";
159                    cout.setf(ios::scientific, ios::floatfield);
160                    cout.precision(15);
161                    for (int j=0; j<numDim; j++)
162                        cout << " " << Nodes->Coordinates[INDEX2(j,i,numDim)];
163                    cout << endl;
164              }              }
165          }          }
166      } else {      } else {
167          fprintf(stdout, "\tNodes: 0D-Nodes 0\n");          cout << "\tNodes: 0D-Nodes 0\n";
168      }      }
169    
170      /* write elements: */      // write elements
171      if (in->Elements != NULL) {      printElementInfo(Elements, "Elements", "Tet4", full);
         int mine = 0, overlap = 0;  
         for (i = 0; i < in->Elements->numElements; i++) {  
             if (in->Elements->Owner[i] == in->MPIInfo->rank)  
                 mine++;  
             else  
                 overlap++;  
         }  
         fprintf(stdout, "\tElements: %s %d (TypeId=%d) owner=%d overlap=%d\n",  
                 in->Elements->ename /*referenceElementSet->referenceElement->Type->Name */ , in->Elements->numElements,  
                 in->Elements->etype /*referenceElementSet->referenceElement->Type->TypeId */ , mine, overlap);  
         NN = in->Elements->numNodes;  
         if (full) {  
             fprintf(stdout, "\t     Id   Tag Owner Color:  Nodes\n");  
             for (i = 0; i < in->Elements->numElements; i++) {  
                 fprintf(stdout, "\t  %5d %5d %5d %5d: ", in->Elements->Id[i], in->Elements->Tag[i],  
                         in->Elements->Owner[i], in->Elements->Color[i]);  
                 for (j = 0; j < NN; j++)  
                     fprintf(stdout, " %5d", in->Nodes->Id[in->Elements->Nodes[INDEX2(j, i, NN)]]);  
                 fprintf(stdout, "\n");  
             }  
         }  
     } else {  
         fprintf(stdout, "\tElements: Tet4 0\n");  
     }  
172    
173      /* write face elements: */      // write face elements
174      if (in->FaceElements != NULL) {      printElementInfo(FaceElements, "Face elements", "Tri3", full);
         int mine = 0, overlap = 0;  
         for (i = 0; i < in->FaceElements->numElements; i++) {  
             if (in->FaceElements->Owner[i] == in->MPIInfo->rank)  
                 mine++;  
             else  
                 overlap++;  
         }  
         fprintf(stdout, "\tFace elements: %s %d (TypeId=%d) owner=%d overlap=%d\n",  
                 in->FaceElements->ename /*referenceElementSet->referenceElement->Type->Name */ ,  
                 in->FaceElements->numElements,  
                 in->FaceElements->etype /*->referenceElementSet->referenceElement->Type->TypeId*/ , mine, overlap);  
         NN = in->FaceElements->numNodes;  
         if (full) {  
             fprintf(stdout, "\t     Id   Tag Owner Color:  Nodes\n");  
             for (i = 0; i < in->FaceElements->numElements; i++) {  
                 fprintf(stdout, "\t  %5d %5d %5d %5d: ", in->FaceElements->Id[i], in->FaceElements->Tag[i],  
                         in->FaceElements->Owner[i], in->FaceElements->Color[i]);  
                 for (j = 0; j < NN; j++)  
                     fprintf(stdout, " %5d", in->Nodes->Id[in->FaceElements->Nodes[INDEX2(j, i, NN)]]);  
                 fprintf(stdout, "\n");  
             }  
         }  
     } else {  
         fprintf(stdout, "\tFace elements: Tri3 0\n");  
     }  
175    
176      /* write points: */      // write points
177      if (in->Points != NULL) {      printElementInfo(Points, "Points", "Point1", full);
         int mine = 0, overlap = 0;  
         for (i = 0; i < in->Points->numElements; i++) {  
             if (in->Points->Owner[i] == in->MPIInfo->rank)  
                 mine++;  
             else  
                 overlap++;  
         }  
         fprintf(stdout, "\tPoints: %s %d (TypeId=%d) owner=%d overlap=%d\n",  
                 in->Points->ename /*->referenceElementSet->referenceElement->Type->Name*/ , in->Points->numElements,  
                 in->Points->etype /*referenceElementSet->referenceElement->Type->TypeId */ , mine, overlap);  
         if (full) {  
             fprintf(stdout, "\t     Id   Tag Owner Color:  Nodes\n");  
             for (i = 0; i < in->Points->numElements; i++)  
             {  
                 fprintf(stdout, "\t  %5d %5d %5d %5d %5d\n", in->Points->Id[i], in->Points->Tag[i],  
                         in->Points->Owner[i], in->Points->Color[i], in->Nodes->Id[in->Points->Nodes[INDEX2(0, i, 1)]]);  
             }  
         }  
     } else {  
         fprintf(stdout, "\tPoints: Point1 0\n");  
     }  
178    
179      /* write tags: */      // write tags
180      if (tag_map) {      if (tagMap.size() > 0) {
181          fprintf(stdout, "\tTags:\n");          cout << "\tTags:\n";
182          while (tag_map) {          TagMap::const_iterator it;
183              fprintf(stdout, "\t  %5d %s\n", tag_map->tag_key, tag_map->name);          for (it=tagMap.begin(); it!=tagMap.end(); it++) {
184              tag_map = tag_map->next;              cout << "\t" << setw(7) << it->second << " " << it->first << endl;
185          }          }
186      }      }
187  }  }

Legend:
Removed from v.6078  
changed lines
  Added in v.6079

  ViewVC Help
Powered by ViewVC 1.1.26