Moved opm/core/vag_format/ -> opm/core/io/vag/

opm/core/io/vag/vag.cpp

@@ -0,0 +1,407 @@
+/*===========================================================================
+//
+// File: vag.cpp
+//
+// Created: 2012-06-08 15:45:53+0200
+//
+// Authors: Knut-Andreas Lie      <Knut-Andreas.Lie@sintef.no>
+//          Halvor M. Nilsen      <HalvorMoll.Nilsen@sintef.no>
+//          Atgeirr F. Rasmussen  <atgeirr@sintef.no>
+//          Xavier Raynaud        <Xavier.Raynaud@sintef.no>
+//          Bård Skaflestad       <Bard.Skaflestad@sintef.no>
+//
+//==========================================================================*/
+
+
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+  Copyright 2012 Statoil ASA.
+
+  This file is part of the Open Porous Media Project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/core/io/vag/vag.hpp>
+#include <opm/core/grid/cornerpoint_grid.h>	 	      
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <stdlib.h>
+#include <cmath>
+#include <cassert>
+#include <set>
+#include <vector>
+#include <map>
+namespace OPM
+{
+    void readPosStruct(std::istream& is,int n,PosStruct& pos_struct){
+	using namespace std;
+	//PosStruct pos_struct;
+	pos_struct.pos.resize(n+1);
+	pos_struct.pos[0]=0;
+	for(int i=0;i< n;++i){
+	    int number;
+	    is >> number ;
+	    //cout <<number << endl;
+	    pos_struct.pos[i+1]=pos_struct.pos[i]+number;
+	    for(int j=0;j< number;++j){
+		int value;
+		is >> value;
+		//	cout << value << " ";
+		pos_struct.value.push_back(value);
+	    }
+	    //cout << endl;
+	}
+	if(int(pos_struct.value.size()) != pos_struct.pos[n]){
+	    cerr << "Failed to read pos structure" << endl;
+	    cerr << "pos_struct.value.size()" << pos_struct.value.size() << endl;
+	    cerr << "pos_struct.pos[n+1]" << pos_struct.pos[n] << endl;
+	}
+    }
+    void writePosStruct(std::ostream& os,PosStruct& pos_struct){
+	using namespace std;
+	//PosStruct pos_struct;
+        if(pos_struct.pos.size()==0){
+            return;
+        }
+        int n=pos_struct.pos.size()-1;        
+	pos_struct.pos.resize(n+1);
+	pos_struct.pos[0]=0;
+	for(int i=0;i< n;++i){
+	    int number=pos_struct.pos[i+1]-pos_struct.pos[i];
+	    os << number << " ";
+	    for(int j=0;j< number;++j){
+		os << pos_struct.value[pos_struct.pos[i]+j] << " ";
+	    }
+	    os << endl;
+	}
+    }
+    void readVagGrid(std::istream& is,OPM::VAG& vag_grid){
+	using namespace std;
+	using namespace OPM;
+	while (!is.eof()) {
+	    string keyword;
+	    is >> keyword;	
+	    //cout << keyword<< endl;
+	    if(keyword == "Number"){
+		string stmp;
+		is >> stmp;
+		if(stmp == "of"){
+		    string entity;
+		    is >> entity;
+		    getline(is,stmp);
+		    int number;
+		    is >> number;
+		    if(entity=="vertices"){
+			vag_grid.number_of_vertices=number;
+		    }else if((entity=="volumes") || (entity=="control")){
+			vag_grid.number_of_volumes=number;
+		    }else if(entity=="faces"){
+			vag_grid.number_of_faces=number;
+		    }else if(entity=="edges"){
+			vag_grid.number_of_edges=number;
+		    }	
+		    cout << "Found Number of: " << entity <<" " << number << endl;		
+		} else {
+		    cerr << "Wrong format: Not of after Number" << endl;
+		    return;
+		}
+	    }else{
+		// read geometry defined by vertices
+		if(keyword=="Vertices"){
+		    int number;
+		    is >> number;	   	    
+		    vag_grid.vertices.resize(3*number);// assume 3d data
+		    readVector(is,vag_grid.vertices);
+		}
+		// here starts the reding of all pos structures
+		else if(keyword=="Volumes->Faces" || keyword=="Volumes->faces"){
+		    //vag_grid.volumes_to_faces=
+		    int number;
+		    is >> number;	
+		    readPosStruct(is,number,vag_grid.volumes_to_faces);
+		    cout << "Volumes->Faces: Number of " << number << endl;
+		}else if(keyword=="Faces->edges" || keyword=="Faces->Edges" ||  keyword=="Faces->Edgess"){
+		    int number;
+		    is >> number;	   
+		    //vag_grid.volumes_to_faces=
+		    readPosStruct(is,number,vag_grid.faces_to_edges);
+		    cout << "Faces->edges: Number of " << number << endl;
+		}else if(keyword=="Faces->Vertices" || keyword=="Faces->vertices"){
+		    int number;
+		    is >> number;	   
+		    //vag_grid.volumes_to_faces=
+		    readPosStruct(is,number,vag_grid.faces_to_vertices);
+		    cout << "Faces->Vertices: Number of " << number << endl;   
+		}else if(keyword=="Volumes->Vertices" || keyword=="Volumes->Verticess"){
+		    int number;
+		    is >> number;	   
+		    //vag_grid.volumes_to_faces=
+		    readPosStruct(is,number,vag_grid.volumes_to_vertices);
+		    cout << "Volumes->Vertices: Number of " << number << endl;   
+		}
+	    
+		//  read simple mappings
+		else if(keyword=="Edge" || keyword=="Edges"){
+		    int number;
+		    is >> number;	   	    
+		    vag_grid.edges.resize(2*number);
+		    readVector(is,vag_grid.edges);
+		    cout << "Edges: Number of " << number << endl;
+		}else if(keyword=="Faces->Volumes" || keyword=="Faces->Control"){
+		    int number;
+		    if(keyword=="Faces->Control"){
+			string vol;
+			is >> vol;
+		    }
+		    is >> number;
+		    vag_grid.faces_to_volumes.resize(2*number);
+		    readVector(is,vag_grid.faces_to_volumes);
+		    cout << "Faces->Volumes: Number of " << number << endl;   
+		}
+		// read material
+		else if(keyword=="Material"){
+		    string snum;
+		    is >> snum;
+		    int number;
+		    is >> number;
+		    cout << "Material number  " << number << endl;
+		    // we read all the rest into doubles		
+		    while(!is.eof()){
+			double value;
+			is >> value;
+			//cout << value << endl;
+			vag_grid.material.push_back(value);
+		    }
+		}else{
+		    //cout << "keyword;
+		}
+		//cout << "Found" << keyword << "Number of " << number << endl;
+	    }
+	}
+    }
+    void vagToUnstructuredGrid(OPM::VAG& vag_grid,UnstructuredGrid& grid){
+	using namespace std;
+	using namespace OPM;
+	cout << "Converting grid" << endl;
+        cout << "Warning:: orignial grid may not be edge confomal" << endl;
+        cout << "          inverse mappings from edges will be wrong" << endl;
+	grid.dimensions=3;
+	grid.number_of_cells=vag_grid.number_of_volumes;
+	grid.number_of_faces=vag_grid.number_of_faces;
+	grid.number_of_nodes=vag_grid.number_of_vertices;
+        
+	// fill face_nodes
+	for(int i=0;i< int(vag_grid.faces_to_vertices.pos.size());++i){
+	    grid.face_nodepos[i] = vag_grid.faces_to_vertices.pos[i];
+	}	    
+	for(int i=0;i< int(vag_grid.faces_to_vertices.value.size());++i){
+	    grid.face_nodes[i] = vag_grid.faces_to_vertices.value[i]-1;
+	}
+	// fill cell_face
+	for(int i=0;i< int(vag_grid.volumes_to_faces.pos.size());++i){
+	    grid.cell_facepos[i] = vag_grid.volumes_to_faces.pos[i];
+	}	    
+	for(int i=0;i< int(vag_grid.volumes_to_faces.value.size());++i){
+	    grid.cell_faces[i] = vag_grid.volumes_to_faces.value[i]-1;
+	}
+	// fill face_cells
+	for(int i=0;i< int(vag_grid.faces_to_volumes.size());++i){
+	    grid.face_cells[i] = vag_grid.faces_to_volumes[i]-1;
+	}
+
+	// fill node_cordinates. This is the only geometry given in the vag
+	for(int i=0;i< int(vag_grid.vertices.size());++i){
+	    grid.node_coordinates[i] = vag_grid.vertices[i];
+	}
+	// informations in edges, faces_to_eges, faces_to_vertices, volume_to_vertices and materials
+	// is not used
+	cout << "Computing geometry" << endl;
+	compute_geometry(&grid);
+	
+    }
+
+    void unstructuredGridToVag(UnstructuredGrid& grid,OPM::VAG& vag_grid){
+        using namespace std;
+        using namespace OPM;
+        cout << "Converting grid" << endl;
+        //	grid.dimensions=3;
+        vag_grid.number_of_volumes=grid.number_of_cells;
+	vag_grid.number_of_faces=grid.number_of_faces;
+	vag_grid.number_of_vertices=grid.number_of_nodes;
+	
+        // resizing vectors
+        vag_grid.vertices.resize(grid.number_of_nodes*3);
+        vag_grid.faces_to_vertices.pos.resize(grid.number_of_faces+1);
+        vag_grid.faces_to_vertices.value.resize(grid.face_nodepos[grid.number_of_faces]);
+        vag_grid.faces_to_volumes.resize(2*grid.number_of_faces);
+        vag_grid.volumes_to_faces.pos.resize(grid.number_of_cells+1);
+        vag_grid.volumes_to_faces.value.resize(grid.cell_facepos[grid.number_of_cells]);//not known
+
+
+        
+        
+        // fill face_nodes
+	for(int i=0;i< int(vag_grid.faces_to_vertices.pos.size());++i){
+	    vag_grid.faces_to_vertices.pos[i] = grid.face_nodepos[i];
+	}
+
+	for(int i=0;i< int(vag_grid.faces_to_vertices.value.size());++i){
+	    vag_grid.faces_to_vertices.value[i] = grid.face_nodes[i] +1;
+	}
+        
+	// fill cell_face      
+	for(int i=0;i< int(vag_grid.volumes_to_faces.pos.size());++i){
+	    vag_grid.volumes_to_faces.pos[i] = grid.cell_facepos[i];
+	}	    
+	for(int i=0;i< int(vag_grid.volumes_to_faces.value.size());++i){
+	    vag_grid.volumes_to_faces.value[i] = grid.cell_faces[i] +1;
+	}
+	// fill face_cells
+	for(int i=0;i< int(vag_grid.faces_to_volumes.size());++i){
+	    vag_grid.faces_to_volumes[i] = grid.face_cells[i] +1;
+	}
+
+	// fill node_cordinates. This is the only geometry given in the vag
+	for(int i=0;i< int(vag_grid.vertices.size());++i){
+	    vag_grid.vertices[i] = grid.node_coordinates[i];
+	}
+
+        
+        // The missing field need to be constructed
+        // gennerate volume to vertice mapping
+        std::vector< std::set<int> > volumes_to_vertices(grid.number_of_cells);
+        for(int i=0;i < grid.number_of_cells; ++i){
+            int nlf=grid.cell_facepos[i+1]-grid.cell_facepos[i];
+            std::set<int> nodes;
+            for(int j=0; j < nlf; ++j){
+                int face = grid.cell_faces[grid.cell_facepos[i]+j];
+                int nlv = grid.face_nodepos[face+1]-grid.face_nodepos[face];
+                for(int k=0; k< nlv; ++k){
+                    int node = grid.face_nodes[grid.face_nodepos[face]+k]+1;
+                    nodes.insert(node);
+                }
+            }
+            volumes_to_vertices[i]=nodes;
+        }
+        
+        // fill volume to vertice map
+        vag_grid.volumes_to_vertices.pos.resize(grid.number_of_cells+1);
+        vag_grid.volumes_to_vertices.value.resize(0);
+        vag_grid.volumes_to_vertices.pos[0]=0;
+        for(int i=0;i < grid.number_of_cells;++i){
+            int nv=volumes_to_vertices[i].size();
+            vag_grid.volumes_to_vertices.pos[i+1]=vag_grid.volumes_to_vertices.pos[i]+nv;
+            std::set<int>::iterator it;
+            for(it=volumes_to_vertices[i].begin();it!=volumes_to_vertices[i].end();++it){
+                vag_grid.volumes_to_vertices.value.push_back(*it);
+            }
+        }
+
+        std::set< std::set<int> > edges;
+	std::vector< std::vector< std::set<int> > > faces_spares;
+	int nfe=0;
+	faces_spares.resize(grid.number_of_faces);
+        for(int i=0;i < grid.number_of_faces;++i){
+            int ne=grid.face_nodepos[i+1]-grid.face_nodepos[i];
+	    nfe=nfe+ne;
+            
+            for(int j=0; j < ne-1;++j){
+                int node1=grid.face_nodes[grid.face_nodepos[i]+j]+1;
+                int node2=grid.face_nodes[grid.face_nodepos[i]+j+1]+1;
+                std::set<int> spair;
+                spair.insert(node1);
+                spair.insert(node2);    
+                edges.insert(spair);
+		faces_spares[i].push_back(spair);     
+            }
+            // add end segment
+            {
+                std::set<int> spair;
+                int node1=grid.face_nodes[grid.face_nodepos[i]+ne-1]+1;
+                int node2=grid.face_nodes[grid.face_nodepos[i]]+1;
+                spair.insert(node1);
+                spair.insert(node2);    
+                edges.insert(spair);
+                faces_spares[i].push_back(spair);     
+            }
+        }
+        
+	// make edge numbering and fill edges
+	std::map<std::set<int>, int> edge_map;
+	std::set< std::set<int> >::iterator it;
+	vag_grid.edges.resize(0);
+	int k=0;
+	for(it=edges.begin(); it!=edges.end();++it){
+	  edge_map.insert(std::pair< std::set<int> , int >(*it,k));
+	  k=k+1;
+	  std::set<int>::iterator sit;
+	  for(sit=(*it).begin();sit!=(*it).end();++sit){
+	    vag_grid.edges.push_back(*sit);
+	  }
+        }
+        // fill face_to_egdes
+	vag_grid.number_of_edges=edges.size();
+	vag_grid.faces_to_edges.pos.resize(vag_grid.number_of_faces+1);
+	for(int i=0;i < grid.number_of_faces;++i){
+            int ne=grid.face_nodepos[i+1]-grid.face_nodepos[i];
+	    vag_grid.faces_to_edges.pos[i+1]=vag_grid.faces_to_edges.pos[i]+ne;	    
+	    for(int j=0;j<int(faces_spares[i].size());++j){
+	      int edge_num=edge_map[faces_spares[i][j]];
+	      vag_grid.faces_to_edges.value.push_back(edge_num+1);
+	    }   
+	}   
+        // vag_grid.edges(0);//not known
+        // vag_grid.faces_to_edges// not known
+        
+        //  material // can not be extracted from the grid
+    }
+    
+    void writeVagFormat(std::ostream& os,OPM::VAG& vag_grid){
+	using namespace std;
+	os << "File in the Vag grid format\n";
+        os << "Number of vertices " ;
+        os << vag_grid.number_of_vertices << endl;;
+        os <<"Number of control volume ";
+        os << vag_grid.number_of_volumes << endl;
+        os <<"Number of faces " ;
+        os << vag_grid.number_of_faces << endl;
+        os <<"Number of edges " ;
+        os << vag_grid.number_of_edges << endl;
+        os <<"Vertices      "   << vag_grid.vertices.size() << endl;
+        writeVector(os, vag_grid.vertices,3);
+        os << "Volumes->faces   " << vag_grid.volumes_to_faces.pos.size()-1 << endl;
+        writePosStruct(os, vag_grid.volumes_to_faces);        
+        os << "Volumes->Vertices   " << vag_grid.volumes_to_vertices.pos.size()-1 << endl;
+        writePosStruct(os, vag_grid.volumes_to_vertices);       
+        os << "Faces->edges   " << vag_grid.faces_to_edges.pos.size()-1 << endl;
+        writePosStruct(os, vag_grid.faces_to_edges);
+        os << "Faces->vertices   " << vag_grid.faces_to_vertices.pos.size()-1 << endl;
+        writePosStruct(os, vag_grid.faces_to_vertices);
+        os << "Faces->Control volumes   " << floor(vag_grid.faces_to_volumes.size()/2) << endl;
+        writeVector(os,vag_grid.faces_to_volumes,2);        
+        os << "Edges   " << floor(vag_grid.edges.size()/2) << endl;
+        writeVector(os,vag_grid.edges,2);
+        /*
+        assert(vag_grid.material.size()%vag_grid.number_of_volumes==0);
+        int lines= floor(vag_grid.material.size()/vag_grid.number_of_volumes);  
+        os << "Material number   " << 1 << endl;
+        writeVector(os,vag_grid.material,lines);
+        */
+
+    }
+
+
+}

opm/core/io/vag/vag.hpp

@@ -0,0 +1,165 @@
+/*===========================================================================
+//
+// File: vag.hpp
+//
+// Created: 2012-06-08 15:46:23+0200
+//
+// Authors: Knut-Andreas Lie      <Knut-Andreas.Lie@sintef.no>
+//          Halvor M. Nilsen      <HalvorMoll.Nilsen@sintef.no>
+//          Atgeirr F. Rasmussen  <atgeirr@sintef.no>
+//          Xavier Raynaud        <Xavier.Raynaud@sintef.no>
+//          Bård Skaflestad       <Bard.Skaflestad@sintef.no>
+//
+//==========================================================================*/
+
+
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+  Copyright 2012 Statoil ASA.
+
+  This file is part of the Open Porous Media Project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_VAG_HPP_HEADER
+#define OPM_VAG_HPP_HEADER
+
+
+
+#include <opm/core/grid.h>
+#include <istream>
+#include <ostream>
+#include <string>
+#include <vector>
+
+namespace OPM
+{
+    /**
+       Struct to hold maping from the natural number less than pos.size()-1 to
+       a set of integers. value(pos(i):pos(i+1)-1) hold the integers corresponding to i.
+       pos(end)-1==value.size();
+    */               
+    struct PosStruct{
+	std::vector<int> pos;
+        std::vector<int>  value;
+    };
+    /**
+      Structure to represent the unstructured vag grid format. The format is only for
+      3D grids.
+     */
+    struct VAG{
+	int number_of_vertices;        
+	int number_of_volumes;
+	int number_of_faces;
+	int number_of_edges;
+        /** Vertices. The coordinates of vertice i is [vetices[3*i:3*i+2]*/
+	std::vector<double> vertices;
+        /** Mapping from volumes to faces */
+	PosStruct volumes_to_faces;
+        /** Mapping from volumes to vertices */
+	PosStruct volumes_to_vertices;
+        /** Mapping from faces to edges */
+	PosStruct faces_to_edges;
+        /** Mapping from faces to vertices */
+	PosStruct faces_to_vertices;
+        /** The edge i is given by the nodes edges[2*i:2*i+1] */
+	std::vector<int> edges;
+        /**  The two neigbours of the face i is faces_to_volumes[2*i:2*i+1] */
+	std::vector<int> faces_to_volumes;
+        /** A vector containing information of each volume. The size is n*number_of_volumes.
+            For each i this is the information:
+            material[n*i] is the volume number and should be transformed to integer
+            material[n*i+1] is a tag and should  be transformed to integer
+            material[n*i+2:n*(i+1)-1] represent propertices.
+        */
+	std::vector<double> material;
+    };
+    /**
+       Function the vag grid format and make a vag_grid struct. This structure
+       is intended to be converted to a grid.
+       \param[in]  is is is stream of the file.
+       \param[out] vag_grid is a reference to a vag_grid struct.
+    */
+    void readVagGrid(std::istream& is,OPM::VAG& vag_grid);
+    /**
+       Function to write vag format.
+       \param[out]  is is is stream of the file.
+       \param[in] vag_grid is a reference to a vag_grid struct.
+
+    */
+    void writeVagFormat(std::ostream& os,OPM::VAG& vag_grid);
+    /**
+       Function to read a vector of some type from a stream.
+       \param[in]  os is is stream of the file.
+       \param[out] vag_grid is a resized and filled vector containing the quantiy read.
+    */   
+    template <typename T>
+    void readVector(std::istream& is,std::vector<T>& vec){
+	using namespace std;
+	for(int i=0;i< int(vec.size());++i){
+	    is >> vec[i];
+	}
+    }
+    /**
+       Function to write a vector of some type from a stream.
+       \param[in]  os is is stream of the file.
+       \param[out] vag_grid is a resized and filled vector containing the quantiy read.
+       \param[in]  n number of doubles on each line.
+    */  
+    template <typename T>
+    void writeVector(std::ostream& os,std::vector<T>& vec,int n){
+        typedef typename std::vector<T>::size_type sz_t;
+
+        const sz_t nn = n;
+
+        for (sz_t i = 0; i < vec.size(); ++i) {
+            os << vec[i] << (((i % nn) == 0) ? '\n' : ' ');
+        }
+
+        if ((vec.size() % nn) != 0) {
+            os << '\n';
+        }
+    }
+
+    /**
+       Read pos struct type mapping from a stream
+       \param[in] is is stream
+       \param[in] n number of lines to read
+       \param[out] pos_struct reference to PosStruct
+     */
+    void readPosStruct(std::istream& is,int n,PosStruct& pos_struct);
+    /**
+       Read pos struct type mapping from a stream
+       \param[in] os is stream to write to
+       \param[in] pos_struct to write
+     */
+    void writePosStruct(std::ostream& os,PosStruct& pos_struct);
+    
+    /**
+       Fill a UnstructuredGrid from a vag_grid.
+       \param[out] vag_grid s is a valid vag_grid struct.
+       \param[in] grid is a grid with have allocated correct size to each pointer.
+     */
+    void vagToUnstructuredGrid(OPM::VAG& vag_grid,UnstructuredGrid& grid);
+    
+    /**
+       Fill a  vag_grid from UnstructuredGrid 
+       \param[out] vag_grid s is a valid vag_grid struct.
+       \param[in] grid is a grid with have allocated correct size to each pointer.
+     */    
+    void unstructuredGridToVag(UnstructuredGrid& grid, OPM::VAG& vag_grid);
+}
+#endif  /* OPM_VAG_HPP_HEADER */	
+