from T2GEORES import geometry as geometry
from T2GEORES import lloyd_relaxation as Field
import numpy as np
import re as re
import subprocess
import datetime
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import shutil
import os
import itertools
import json
import shapefile
import pylab as plb
import math
import sys
from scipy.spatial import ConvexHull
from scipy.interpolate import griddata
from scipy.spatial.distance import cdist
from scipy.spatial import cKDTree
import pandas as pd
from matplotlib.ticker import (MultipleLocator, FormatStrFormatter,AutoMinorLocator)
import pyvista as pv
import vtk
import sqlite3
import geopandas as gpd
import string
[docs]class py2amesh:
"""It creates a mesh based on well positions and irregular blocks
The main characteristics are:
-It generates the mesh based on defined boundaries
-The inner section is called wellfield and can contain elements with hexagonal or square shape
-The wellfield can be delimited by a shapefile or by fixed squared boundaries.
-It allows to generates elements along a line to represent structures.
-The voronoi elements are generated by AMESH.
-The elements name cannot contain two consecutive zeros.
-The well blocks are the first on the list
-It creates the input files to work on Steinar (RockEditor)
-It can export a defined layer on a shapefile format
-It can plot a selected layer (It is recommended to use the function plot_voronoi() to plot)
-Two json file are generated with the correlative block for each well, from which can be track during all the modelling steps.
Parameters
----------
filename : str
File name with well feedzone location
filepath : str
Path of input files
Xmin : float
Minimun X coordinates for the grid
Xmax : float
Maximun X coordinates for the grid
Ymin : float
Minimun Y coordinates for the grid
Ymax : float
Maximun Y coordinates for the grid
toler : float
AMESH parameter
layers : dictionary
Name (correlative) and thickness of every layer on the model, keyword on input_dictionary
layer_to_plot : int
In case it is specified a voronoi plot will be performed
x_space : float
Horizontal distance between elements for the outerfield
y_space : float
Vertical distance between elements for the outerfield
radius_criteria: float
Minimun distance between well location and a regular element
x_from_boarder: float
Horizontal distance from the first element to the east border
y_from_boarder: float
Vertical distance from the first element to the south border
x_gap_min: float
Minimun X coordinates on the grid for the well field
x_gap_max: float
Maximun X coordinates on the grid for the well field
x_gap_space: float
Horizontal distance between elements for the farfield
y_gap_min: float
Minimun Y coordinates on the grid for the well field o
y_gap_max: float
Maximun X coordinates on the grid for the well field
y_gap_space: float
Vertical distance between elements for the farfield
plot_names: bool
If true it plots the name of the blocks from the selected layer to plot
plot_centers: bool
If true it plots the centers of the blocks from the selected layer to plot
z0_level: float
Reference level (elevation) for all the grid, keyword on input_dictionary
mesh_creation: bool
If true the mesh is created
plot_layer: bool
If true it plots the selected layer
to_steinar: bool
If true it creates the input files for steinar
to_GIS: bool
If true it generates a shapefile of the selected layer
plot_all_GIS: bool
If true it generates a shapefile of all layers
from_leapfrog: bool
lee archivos leapfrong ../mesh/from_leapfrog/LF_geometry.dat y ../mesh/from_leapfrog/LF_t2.dat, sin embargo se pierde la simbologia usada en leapfrog y no te utiliza la malla regular ni los pozos. Solamente se crea la malla utilizando amesh
line_file: str
It defines the path and name of a line that can represented a fault or other structure on the mesh, The input file must contain the header: ID,X,Y on csv format. ID referes to the same
structure, thus, more than one structure can be defined on a single file.
fault_distance: float
In case a line_file is define, some paralels elements will be created at a defined distance
with_polygon: bool
If true a shapefile will be read to define the wellfield.
polygon_shape: str
The shapefile deines the wellfield boundaries. The shape must not contain any cavity
set_inac_from_poly: bool
If true all the elements on the outside of the shapefile are defined as inactive
set_inac_from_inner:bool
If true all the elements on the outerfield are defined as inactive
rotate: bool
If true it rotates the mesh a defined angle
angle: float
Angle in degrees
inner_mesh_type: string
Type of mesh on the inner part of the mesh, it could be 'honeycomb' or 'regular'
Returns
-------
file
eleme: list of blocks from the grid
file
conne : list of connections on the grid
shapefile
mesh_{field}_layer_{layer} : shapefile of a defined (or all) layer including rock distribution
plot
Voronoi plot (in case is specified)
Attention
---------
A copy of AMESH must be on the path or directory
"""
def __init__(self,filename,filepath,Xmin,Xmax,Ymin,Ymax,\
toler,layers,layer_to_plot,x_space,y_space,radius_criteria,\
x_from_boarder,y_from_boarder,\
x_gap_min,x_gap_max,x_gap_space,y_gap_min,y_gap_max,y_gap_space,\
plot_names,plot_centers,z0_level,plot_all_GIS,from_leapfrog,line_file,fault_distance,with_polygon,polygon_shape,set_inac_from_poly,set_inac_from_inner,rotate,angle,inner_mesh_type,\
distance_points,fault_rows,relaxation_times,points_around_well,distance_points_around_well, outer_polygon):
self.filename=filename
self.filepath=filepath
self.layers=layers
self.number_of_layer=len(layers)
self.Xmin=Xmin
self.Xmax=Xmax
self.Ymin=Ymin
self.Ymax=Ymax
self.z0_level=z0_level
self.layer_to_plot=layer_to_plot
self.radius_criteria=radius_criteria
self.x_space=x_space
self.y_space=y_space
self.x_from_boarder=x_from_boarder
self.y_from_boarder=y_from_boarder
self.z=0
self.delf_rock="101" #Layer index
self.filename_out="in"
self.toler=toler
self.x_gap_min=x_gap_min
self.x_gap_max=x_gap_max
self.x_gap_space=x_gap_space
self.y_gap_min=y_gap_min
self.y_gap_max=y_gap_max
self.y_gap_space=y_gap_space
self.plot_names=plot_names
self.plot_centers=plot_centers
self.plot_all_GIS=plot_all_GIS
self.from_leapfrog=from_leapfrog
self.line_file=line_file
self.fault_distance=fault_distance
self.with_polygon=with_polygon
self.set_inac_from_poly=set_inac_from_poly
self.set_inac_from_inner=set_inac_from_inner
self.rotate=rotate
self.angle=angle
self.inner_mesh_type=inner_mesh_type
self.polygon_shape=polygon_shape
self.distance_points=distance_points
self.fault_rows=fault_rows
self.relaxation_times=relaxation_times
self.points_around_well=points_around_well
self.distance_points_around_well=distance_points_around_well
self.outer_polygon = outer_polygon
"""
shape = shapefile.Reader(polygon_shape)
#first feature of the shapefile
feature = shape.shapeRecords()[0]
points = feature.shape.__geo_interface__
self.polygon=[]
for n in points:
for v in points[n]:
if n=='coordinates':
self.polygon.append([v[0],v[1]]) # (GeoJSON format)
"""
inner_borders=gpd.read_file(polygon_shape)
polygon=[]
for line in inner_borders.iterrows():
pointList = line[1].geometry.exterior.coords.xy
for point in zip(pointList[0],pointList[1]):
polygon.append([point[0],point[1]])
self.polygon = polygon[::-1][0:-1]
#Read border to clip write into in files
borders=gpd.read_file(outer_polygon)
border_points=[]
for line in borders.iterrows():
pointList = line[1].geometry.exterior.coords.xy
for point in zip(pointList[0],pointList[1]):
border_points.append([point[0],point[1]])
self.polygon_external=border_points[::-1][0:-1]
self.color_dict = {1:[['AA','AB','AC','AD','AE','AF','AG'],'ROCK1','red'],\
2:[['BA','BB','BC','BD','BE','BF','BG'],'ROCK2','white'],\
3:[['CA','CB','CC','CD','CE','CF','CG'],'ROCK3','yellow'],\
4:[['DA','DB','DC','DD','DE','DF','DG'],'ROCK4','blue'],\
5:[['EA','EB','EC','ED','EE','EF','EG'],'ROCK5','green'],\
6:[['FA','FB','FC','FD','FE','FF','FG'],'ROCK6','purple'],\
7:[['GA','GB','GC','GD','GE','GF','GG'],'ROCK7','#ff69b4'],\
8:[['HA','HB','HC','HD','HE','HF','HG'],'ROCK8','darkorange'],\
9:[['IA','IB','IC','ID','IE','IF','IG'],'ROCK9','cyan'],\
10:[['JA','JB','JC','JD','JE','JF','JG'],'ROK10','magenta'],\
11:[['KA','KB','KC','KD','KE','KF','KG'],'ROK11','#faebd7'],\
12:[['LA','LB','LC','LD','LE','LF','LG'],'ROK12','#2e8b57'],\
13:[['MA','MB','MC','MD','ME','MF','MG'],'ROK13','#eeefff'],\
14:[['NA','NB','NC','ND','NE','NF','NG'],'ROK14','#da70d6'],\
15:[['OA','OB','OC','OD','OE','OF','OG'],'ROK15','#ff7f50'],\
16:[['PA','PB','PC','PD','PE','PF','PG'],'ROK16','#cd853f'],\
17:[['QA','QB','QC','QD','QE','QF','QG'],'ROK17','#bc8f8f'],\
18:[['RA','RB','RC','RD','RE','RF','RG'],'ROK18','#5f9ea0'],\
19:[['SA','SB','SC','SD','SE','SF','SG'],'ROK19','#daa520'],
20:[['TA','TB','SC','SD','SE','SF','SG'],'ROK20','#daa520'],
21:[['UA','UB','UC','UD','UE','UF','UG'],'ROK21','#daa520'],
22:[['VA','VB','SC','VD','VE','VF','VG'],'ROK22','#daa520'],
23:[['WA','WB','SC','WD','WE','WF','WG'],'ROK23','#daa520'],
24:[['XA','XB','SC','XD','XE','XF','XG'],'ROK19','#daa520'],
25:[['YA','YB','SC','YD','YE','YF','YG'],'ROK19','#daa520'],
26:[['ZA','ZB','SC','ZD','ZE','ZF','ZG'],'ROK19','#daa520']}
self.rock_dict={}
prof_cont=0
for jk in range(1,len(layers)+1):
if jk==1:
prof_cont=layers[jk-1]*0.5
z_real=z0_level-prof_cont
elif jk>1:
prof_cont=prof_cont+layers[jk-1]*0.5+layers[jk-2]*0.5
z_real=z0_level-prof_cont
self.rock_dict[jk]=[self.color_dict[jk][0][0],self.color_dict[jk][1],\
self.color_dict[jk][2],self.color_dict[jk][0][0],z_real,self.layers[jk-1]]
[docs] def regular_mesh(self):
"""Genera malla regular en en toda la extension de la region definida por Xmin,Xmax,Ymin y Ymax
"""
x_regular=range(self.Xmin+self.x_from_boarder,self.Xmax+self.x_space-self.x_from_boarder,self.x_space)
y_regular=range(self.Ymin+self.y_from_boarder,self.Ymax+self.y_space-self.y_from_boarder,self.y_space)
x_regular_small=range(self.x_gap_min,self.x_gap_max+self.x_gap_space,self.x_gap_space)
y_regular_small=range(self.y_gap_min,self.y_gap_max+self.y_gap_space,self.y_gap_space)
self.mesh_array=[]
for nx in x_regular:
for ny in y_regular:
if ((nx<self.x_gap_min) or (nx>self.x_gap_max)) or ((ny<self.y_gap_min) or (ny>self.y_gap_max)):
self.mesh_array.append([nx,ny])
#Small polygon area must be here
for nxx in x_regular_small:
cnt=0
for nyy in y_regular_small:
if [nxx,nyy] not in self.mesh_array:
if self.inner_mesh_type=='honeycomb':
if cnt%2==0:
self.mesh_array.append([nxx,nyy])
else:
self.mesh_array.append([nxx+self.x_gap_space/2,nyy])
elif self.inner_mesh_type=='regular':
self.mesh_array.append([nxx,nyy])
cnt+=1
if self.rotate:
angle=self.angle
for pair in range(len(self.mesh_array)):
x1=self.mesh_array[pair][0]-self.Xmin
y1=self.mesh_array[pair][1]-self.Ymin
self.mesh_array[pair][0]=x1*math.cos(math.pi*angle/180)-y1*math.sin(math.pi*angle/180)+self.Xmin
self.mesh_array[pair][1]=x1*math.sin(math.pi*angle/180)+y1*math.cos(math.pi*angle/180)+self.Ymin
return np.array(self.mesh_array)
[docs] def check_in_out(self,position,point,source):
"""Verifica si un punto de la malla del campo cercano esta dentro o fuera del poligo definido por el shapefile de entrada o del campo cercano
"""
if position=='internal':
polygon=self.polygon
elif position=='external':
polygon=self.polygon_external
boolean=False
if source=='shapefile':
cnt=0
for n in range(len(polygon)):
if n+1!=len(polygon):
if polygon[n][0]==polygon[n+1][0]:
m = np.inf
b = np.inf
x = polygon[n][0]
elif polygon[n][1] == polygon[n+1][1]:
m = 0
b = polygon[n][1]
x = point[0]
else:
m,b=plb.polyfit([polygon[n][0],polygon[n+1][0]],[polygon[n][1],polygon[n+1][1]],1)
val_range=[polygon[n][1],polygon[n+1][1]]
elif n+1==len(polygon):
if polygon[-1][0] == polygon[0][0]:
m = np.inf
b = np.inf
x = polygon[-1][0]
elif polygon[-1][1] == polygon[0][1]:
m = 0
b = polygon[-1][1]
x = point[0]
else:
m,b=plb.polyfit([polygon[-1][0],polygon[0][0]],[polygon[-1][1],polygon[0][1]],1)
val_range=[polygon[-1][1],polygon[0][1]]
if m != np.inf and m != 0:
x = (point[1]-b)/m
if point[0]<x and min(val_range)<point[1] and point[1]<max(val_range):
cnt+=1
if cnt==1:
boolean=True
elif source=='inner':
Xarray_inner=np.array([self.x_gap_min,self.x_gap_max+self.x_gap_space])
Yarray_inner=np.array([self.y_gap_min,self.y_gap_max+self.y_gap_space])
if Xarray_inner[0]<point[0] and Xarray_inner[1]>point[0] and Yarray_inner[0]<point[1] and Yarray_inner[1]>point[1]:
boolean=True
return boolean
[docs] def reg_pol_mesh(self):
"""Crea malla regular cuando existe un poligono de entrada
"""
#x_regular=np.arange(self.Xmin+self.x_from_boarder,self.Xmax+self.x_space-self.x_from_boarder,self.x_space)
#y_regular=np.arange(self.Ymin+self.y_from_boarder,self.Ymax+self.y_space-self.y_from_boarder,self.y_space)
nx=40 #number of elements in one direction
n_times=4
ny=int((self.Ymax+2*n_times*self.x_from_boarder-self.Ymin)*nx/(self.Xmax-self.Xmin+n_times*2*self.x_from_boarder))
x_regular=np.linspace(self.Xmin-n_times*self.x_from_boarder,self.Xmax+n_times*self.x_from_boarder,nx,endpoint=True) #generation of regular grid on X
y_regular=np.linspace(self.Ymin-n_times*self.y_from_boarder,self.Ymax+n_times*self.x_from_boarder,ny,endpoint=True) #generation of regular grid on Y
self.mesh_array=[]
for nx in x_regular:
for ny in y_regular:
self.mesh_array.append([nx,ny])
if self.rotate:
angle=self.angle
for pair in range(len(self.mesh_array)):
x1=self.mesh_array[pair][0]-(self.Xmin-n_times*self.x_from_boarder)
y1=self.mesh_array[pair][1]-(self.Ymin-n_times*self.x_from_boarder)
self.mesh_array[pair][0]=x1*math.cos(math.pi*angle/180)-y1*math.sin(math.pi*angle/180)+self.Xmin-n_times*self.x_from_boarder
self.mesh_array[pair][1]=x1*math.sin(math.pi*angle/180)+y1*math.cos(math.pi*angle/180)+self.Ymin-n_times*self.x_from_boarder
x_pol=[]
y_pol=[]
for n in range(len(self.polygon)):
x_pol.append(int(self.polygon[n][0]))
y_pol.append(int(self.polygon[n][1]))
x_pol.append(int(self.polygon[0][0]))
y_pol.append(int(self.polygon[0][1]))
x_gap_min=self.x_gap_min#min(x_pol)
x_gap_max=self.x_gap_max#max(x_pol)
y_gap_min=self.y_gap_min#min(y_pol)
y_gap_max=self.y_gap_max#max(y_pol)
small_mesh=[]
#x_regular_small=np.arange(x_gap_min,x_gap_max+self.x_gap_space,self.x_gap_space)
#y_regular_small=np.arange(y_gap_min,y_gap_max+self.y_gap_space,self.y_gap_space)
x_regular_small=np.arange(x_gap_min,x_gap_max,self.x_gap_space)
y_regular_small=np.arange(y_gap_min,y_gap_max,self.y_gap_space)
for nxx in x_regular_small:
cnt=0
for nyy in y_regular_small:
if [nxx,nyy] not in small_mesh:
if self.inner_mesh_type=='honeycomb':
if cnt%2==0:
small_mesh.append([nxx,nyy])
else:
small_mesh.append([nxx+self.x_gap_space/2,nyy])
elif self.inner_mesh_type=='regular':
small_mesh.append([nxx,nyy])
cnt+=1
if self.rotate:
angle=self.angle
for pair in range(len(small_mesh)):
x1=small_mesh[pair][0]-self.x_gap_min
y1=small_mesh[pair][1]-self.y_gap_min
small_mesh[pair][0]=x1*math.cos(math.pi*angle/180)-y1*math.sin(math.pi*angle/180)+self.x_gap_min
small_mesh[pair][1]=x1*math.sin(math.pi*angle/180)+y1*math.cos(math.pi*angle/180)+self.y_gap_min
to_delete=[]
for v in range(len(self.mesh_array)):
point=[self.mesh_array[v][0],self.mesh_array[v][1]]
check=self.check_in_out('internal',point,source='shapefile')
if check:
to_delete.append(v)
self.mesh_array=np.delete(self.mesh_array, to_delete, 0)
to_delete=[]
for v in range(len(small_mesh)):
point=[small_mesh[v][0],small_mesh[v][1]]
check=self.check_in_out('internal',point,source='shapefile')
if not check:
to_delete.append(v)
small_mesh=np.delete(small_mesh, to_delete, 0)
mesh_pol=[]
for vk in range(len(self.mesh_array)):
mesh_pol.append([self.mesh_array[vk][0],self.mesh_array[vk][1]])
for vk in range(len(small_mesh)):
mesh_pol.append([small_mesh[vk][0],small_mesh[vk][1]])
return np.array(mesh_pol)
[docs] def radius_select(self,x0,y0,xr,yr,type_i='mesh'):
"""Verifica si dos puntos estan mas cerca que el criterio seleccionado
"""
r=((x0-xr)**2+(y0-yr)**2)**0.5
if type_i=='mesh':
cx=self.radius_criteria
elif type_i=='well':
cx=10*0.4*2**0.5
elif type_i=='fault':
cx=15
if r<cx:
boolean=1
else:
boolean=0
return boolean
[docs] def from_leapfrog_mesh(self):
"""Extrae puntos mas extremos y la posicion de los elementos de un set de datos provinientes de leapfrog, sin embargo no considera los elementos asociados a la roca ATM 0
"""
geometry_file="../mesh/from_leapfrog/LF_geometry.dat"
leapfrog_t2_file="../mesh/from_leapfrog/LF_t2.dat"
#Creates a dictionary using the layers
printlayer=False
layer_min=[]
layers={}
with open(geometry_file,'r') as f:
for line in f.readlines():
if line.rstrip()=='LAYERS':
printlayer=True
continue
elif line.rstrip()=='SURFA' or line.rstrip()=='':
printlayer=False
if printlayer:
layer=line.rstrip()[0:2]
if layer==' 0':
layer_min.append(line.rstrip()[2:13])
layer_middle=line.rstrip()[13:23]
layer_max=line.rstrip()[13:23]
continue
else:
layer_max=layer_min[-1]
layer_min.append(line.rstrip()[2:13])
layer_middle=line.rstrip()[13:23]
layers[int(layer)]=[float(layer_max),float(layer_middle),float(layer_min[-1])]
max_layer=max(layers.keys())
xc=[]
yc=[]
self.LP_mesh=[]
printeleme=False
#Takes the elements at the selected leyer
with open(leapfrog_t2_file,'r') as f:
for line in f.readlines():
if line.rstrip()=='ELEME':
printeleme=True
continue
elif line.rstrip()=='CONNE' or line.rstrip()=='':
printeleme=False
if printeleme and line.rstrip()[0:5]!="ATM 0" and int(line.rstrip()[3:5])==max_layer:
xc=float(line.rstrip()[51:60])
yc=float(line.rstrip()[60:70])
self.LP_mesh.append([xc,yc])
#Creates a dictionary of the VERTICES
printvertices=False
self.x_min=1E20
self.x_max=0
self.y_min=1E20
self.y_max=0
with open(geometry_file,'r') as f:
for line in f.readlines():
#It will read between the keywords GRID and CONNECTIONS
if line.rstrip()=='VERTICES':
printvertices=True
continue
elif line.rstrip()=='GRID' or line.rstrip()=="":
printvertices=False
if printvertices:
vertice_x=float(line.rstrip()[4:13])
vertice_y=float(line.rstrip()[13:23])
if vertice_x>self.x_max:
self.x_max=vertice_x
if vertice_y>self.y_max:
self.y_max=vertice_y
if vertice_x<self.x_min:
self.x_min=vertice_x
if vertice_y<self.y_min:
self.y_min=vertice_y
return self.x_max,self.x_min, self.y_min,self.y_max, self.LP_mesh
[docs] def data(self):
"""Define los puntos que ingresaran al archivo de entrada para amesh. Adicionalmente, en caso de definir un linea en el archivo de entrada <i><lines_data/i> se procedera a ingresar estos puntos y crear puntos paralelos en ambos extremos de la linea
"""
self.raw_data=np.genfromtxt(self.filepath+self.filename,dtype={'names':('ID','MD','X','Y','Z','TYPE'),'formats':('<U7','f4','f4','f4','f4','<U10')},delimiter=',',skip_header=True)
self.IDXY={}
if not self.from_leapfrog:
if self.with_polygon:
regular_mesh=self.reg_pol_mesh()
else:
regular_mesh=self.regular_mesh()
else:
x,x1,y1,y2,regular_mesh=self.from_leapfrog_mesh()
for n in range(len(self.raw_data['ID'])):
#Store the data from wells
self.IDXY["%s"%(str(self.raw_data['ID'][n]))]=[self.raw_data['X'][n],self.raw_data['Y'][n],self.raw_data['TYPE'][n]]
#self.IDXY["%s"%(str(self.raw_data['ID'][n])).split("'")[1]]=[self.raw_data['X'][n],self.raw_data['Y'][n],self.raw_data['TYPE'][n]]
to_delete=[]
x0=self.raw_data['X'][n]
y0=self.raw_data['Y'][n]
#Delete the regular points close to the wells
for ngrid in range(len(regular_mesh)):
if abs(x0-regular_mesh[ngrid][0])<self.radius_criteria or abs(y0-regular_mesh[ngrid][1])<self.radius_criteria:
boolean=self.radius_select(x0,y0,regular_mesh[ngrid][0],regular_mesh[ngrid][1])
if boolean==1:
to_delete.append(ngrid)
regular_mesh=np.delete(regular_mesh, to_delete, 0)
#Mesh around the wells
if self.points_around_well > 0:
well_mesh=[]
d=self.distance_points_around_well
dw=np.linspace(-d*0.4,d*0.4,num=self.points_around_well)
dw=dw[dw!=0]
for n in range(len(self.raw_data['ID'])):
x=self.raw_data['X'][n]
y=self.raw_data['Y'][n]
for xr in dw:
for yr in dw:
xi=x+xr
yi=y+yr
#Rotating mesh around the wells
if self.rotate:
angle=self.angle
x1=xi-(x)
y1=yi-(y)
xii=x1*math.cos(math.pi*angle/180)-y1*math.sin(math.pi*angle/180)+(x)
yii=x1*math.sin(math.pi*angle/180)+y1*math.cos(math.pi*angle/180)+(y)
well_mesh.append([xii,yii])
else:
well_mesh.append([xi,yi])
well_mesh=np.array(well_mesh)
#Deleting elements on the mesh around the well to close to other regular elements of the mesh
for point in regular_mesh:
to_delete=[]
for i,well_point in enumerate(well_mesh):
boolean=self.radius_select(well_point[0],well_point[1],point[0],point[1],type_i='well')
if boolean==1:
to_delete.append(i)
well_mesh=np.delete(well_mesh, to_delete, 0)
#Finally appending the mesh
regular_mesh=np.append(regular_mesh,well_mesh,axis=0)
cnt=0
if self.line_file:
self.lines_data=np.loadtxt(self.filepath+self.line_file,dtype={'names':('ID','X','Y'),'formats':('S10','f4','f4')},delimiter=',',skiprows=1)
lines_dict={}
for v in self.lines_data:
key=v[0].decode('UTF-8')
if key not in list(lines_dict.keys()):
lines_dict[key]=[[float(v[1]),float(v[2])]]
else:
lines_dict[key].append([float(v[1]),float(v[2])])
d=self.fault_distance #Fine more of this for better definition a long the faults
#ds=np.arange(d,d*2.5,step=d/2)
#df=self.x_gap_space-d/4 #step=d/2
ds=[]
init_d=d/2
for j in range(self.fault_rows):
d_val=init_d*2**(j)
ds.append(d_val)
ds=np.array(ds)
#Iter over every structure
for n in lines_dict:
x_lines=[]
y_lines=[]
for v in lines_dict[n]:
x_lines.append(v[0])
y_lines.append(v[1])
x_lines_np=np.array([]) #Initialize X position of structures
y_lines_np=np.array([]) #Initialize Y position of structures
for nv in range(len(x_lines)):
if nv+1<len(x_lines):
#Generates a line with the usual two points of the fault
line=[(x_lines[nv],y_lines[nv]),(x_lines[nv+1],y_lines[nv+1])]
#Calculates the distance between the points
distDiag = cdist(line,line,'euclidean').diagonal(1)
#Estimastes the number of new points
d_btw_points=self.distance_points # 10
num_steps=int(distDiag//d_btw_points)
#Generates the new x positions to evaluate
x_lines_np_local=np.linspace(x_lines[nv],x_lines[nv+1],num=num_steps)
y_lines_np_local=np.linspace(y_lines[nv],y_lines[nv+1],num=num_steps)
#Generates a double quantity of new x positions
x_lines_np_local2=[]
y_lines_np_local2=[]
for h, val in enumerate(x_lines_np_local):
if h+1<len(x_lines_np_local):
x_val=(val+x_lines_np_local[h+1])/2
y_val=(y_lines_np_local[h]+y_lines_np_local[h+1])/2
x_lines_np_local2.append(x_val)
y_lines_np_local2.append(y_val)
x_lines_np_local2=np.array(x_lines_np_local2)
y_lines_np_local2=np.array(y_lines_np_local2)
x_lines_np=np.append(x_lines_np,x_lines_np_local)
y_lines_np=np.append(y_lines_np,y_lines_np_local)
#calculates line parameters
m,b=plb.polyfit(x_lines_np_local,y_lines_np_local,1)
#perpendicular splope
p_s=(-1/m)
for i, d in enumerate(ds):
c=d*(1+m**2)**0.5+b
if b<0:
c=b-d*(1+m**2)**0.5
if i%2!=0:
b2p=[]
for point in zip(x_lines_np_local,y_lines_np_local):
b2p.append(point[1]-point[0]*p_s)
if b>0:
x_u=[(c-b2)/(p_s-m) for x,b2 in zip(x_lines_np_local,b2p)]
x_d=[(2*b-c-b2)/(p_s-m) for x,b2 in zip(x_lines_np_local,b2p)]
else:
x_u=[(2*b-c-b2)/(p_s-m) for x,b2 in zip(x_lines_np_local,b2p)]
x_d=[(c-b2)/(p_s-m) for x,b2 in zip(x_lines_np_local,b2p)]
y_u=[p_s*xu+b2 for xu,b2 in zip(x_u,b2p)]
y_d=[p_s*xd+b2 for xd,b2 in zip(x_d,b2p)]
x_lines_np=np.append(x_lines_np,x_u)
x_lines_np=np.append(x_lines_np,x_d)
y_lines_np=np.append(y_lines_np,y_u)
y_lines_np=np.append(y_lines_np,y_d)
else:
x_u=[]
x_d=[]
y_u=[]
y_d=[]
for i2, x in enumerate(x_lines_np_local2):
b2p=y_lines_np_local2[i2]-x*p_s
if b>0:
x_uu=(c-b2p)/(p_s-m)
x_dd=(2*b-c-b2p)/(p_s-m)
else:
x_uu=(2*b-c-b2p)/(p_s-m)
x_dd=(c-b2p)/(p_s-m)
x_u.append(x_uu)
x_d.append(x_dd)
y_uu=p_s*x_uu+b2p
y_u.append(y_uu)
y_dd=p_s*x_dd+b2p
y_d.append(y_dd)
x_lines_np=np.append(x_lines_np,x_u)
x_lines_np=np.append(x_lines_np,x_d)
y_lines_np=np.append(y_lines_np,y_u)
y_lines_np=np.append(y_lines_np,y_d)
for n2 in range(len(x_lines_np)):
to_delete=[]
x0=x_lines_np[n2]
y0=y_lines_np[n2]
for ngrid in range(len(regular_mesh)):
boolean=self.radius_select(x0,y0,regular_mesh[ngrid][0],regular_mesh[ngrid][1],type_i='fault')
if boolean==1:
to_delete.append(ngrid)
regular_mesh=np.delete(regular_mesh, to_delete, 0)
for n3 in range(len(x_lines_np)):
self.IDXY["RG%s"%(n3+cnt)]=[x_lines_np[n3],y_lines_np[n3],'REGUL']
cnt+=n3
#Store the data from the regular mesh on IDXY
for n in range(len(regular_mesh)):
self.IDXY["RG%s"%(n+cnt)]=[regular_mesh[n][0],regular_mesh[n][1],'REGUL']
return {'raw_data':self.raw_data,'IDXY':self.IDXY}
[docs] def well_blk_assign(self):
"""Asigna el nombre a cada bloque, priorizando los pozos, es decir, estos llevan los correlativos mas bajos. Por ultimo almacena dos archivo json para registro de los bloques asignados a los pozos
"""
self.well_names=[]
self.blocks_PT={}
self.well_blocks_names={}
self.wells_correlative={}
data_dict=self.data()
for n in range(self.number_of_layer):
cnt=110
layer_num=0
for x in data_dict['IDXY']:
if cnt%100==0:
blocknumber=cnt+10
cnt=blocknumber
else:
blocknumber=cnt
if cnt%1010==0:
cnt=110
blocknumber=cnt
layer_num+=1
cnt+=1
string_ele=self.color_dict[n+1][0][(layer_num)]
blockname=string_ele+str(blocknumber)
self.well_names.append(blockname)
self.well_blocks_names["%s"%(blockname)]=[x,\
data_dict['IDXY'][x][0],\
data_dict['IDXY'][x][1],\
data_dict['IDXY'][x][2],\
n+1,self.rock_dict[n+1][4],\
self.rock_dict[n+1][2],self.rock_dict[n+1][5]]
#print(data_dict['IDXY'][x][2]!='REGUL')
#if data_dict['IDXY'][x][2]=='prod' or data_dict['IDXY'][x][2]=='rein':
if data_dict['IDXY'][x][2]!='REGUL':
#str(data_dict['IDXY'][x][2]).split("'")[1],\
self.blocks_PT["%s"%(blockname)]=[x,\
str(data_dict['IDXY'][x][0]),\
str(data_dict['IDXY'][x][1]),\
str(data_dict['IDXY'][x][2]),\
str(n+1),str(self.rock_dict[n+1][4]),\
self.rock_dict[n+1][2],str(self.rock_dict[n+1][5])]
try:
self.wells_correlative[x]=["%s"%(blockname[1:])]
except KeyError:
pass
json.dump(self.blocks_PT, open("../mesh/well_dict.txt",'w'),sort_keys=True, indent=4)
json.dump(self.wells_correlative, open("../mesh/wells_correlative.txt",'w'),sort_keys=True, indent=4)
return {'well_names':self.well_names, 'well_blocks_names':self.well_blocks_names}
[docs] def to_translate(self,points):
#Converts data points to new coordinate system
n_times=4
if self.rotate:
angle=self.angle
for n, point in enumerate(points):
if point[1]!=(self.Ymin-n_times*self.x_from_boarder) and point[0]!=(self.Xmin-n_times*self.x_from_boarder):
line=[(self.Xmin-n_times*self.x_from_boarder,self.Ymin-n_times*self.x_from_boarder),(point[0],point[1])]
r = cdist(line,line,'euclidean').diagonal(1)
alpha = np.arctan((point[1]-self.Ymin+n_times*self.x_from_boarder)/(point[0]-self.Xmin+n_times*self.x_from_boarder))
if alpha<0:
alpha=np.pi+alpha
betha = alpha - np.deg2rad(angle)
points[n][0]=r*np.cos(betha)
points[n][1]=r*np.sin(betha)
else:
points[n][0]=0
points[n][1]=0
return points
[docs] def de_translate(self,points):
#Converts data points to new coordinate system
n_times=4
if self.rotate:
angle=self.angle
for n, point in enumerate(points):
if point[1]!=0 and point[0]!=0:
betha = np.arctan(point[1]/point[0])
alpha=betha+np.deg2rad(angle)
line=[(0,0),(point[0],point[1])]
r = cdist(line,line,'euclidean').diagonal(1)
points[n][0]=r*np.cos(alpha)+(self.Xmin-n_times*self.x_from_boarder)
points[n][1]=r*np.sin(alpha)+(self.Ymin-n_times*self.x_from_boarder)
else:
points[n][0]=self.Xmin
points[n][1]=self.Ymin
return points
[docs] def well_block_assign_from_relaxation(self,layer,i):
letters=string.ascii_uppercase
name=layer+letters[int(i/810)]+str(int(1+(i-810*int(i/810))/90))+str(10+i-90*int(i/90))
return name
[docs] def relaxation(self):
data={'blocks':{},'borders':[]}
#Extracts the points from the regular mesh, it considers the rotation already and internal region of shapefile
data_dict=self.data()['IDXY']
points=[]
for key in data_dict:
points.append([data_dict[key][0],data_dict[key][1]])
points=np.array(points)
points=np.array(points)
fig, ax = plt.subplots(1, 1, figsize=(20,20))
ax.plot(points[:,0],points[:,1],'ok', linestyle='None',ms=1)
for k, dot in enumerate(self.polygon_external):
if k+1<len(self.polygon_external):
ax.plot([dot[0],self.polygon_external[k+1][0]],[dot[1],self.polygon_external[k+1][1]],'og')
ax.set_aspect('equal')
plt.savefig('points_original.png',dpi=600)
plt.show()
points=self.to_translate(points)
points=np.array(points)
fig, ax = plt.subplots(1, 1, figsize=(20,20))
ax.plot(points[:,0],points[:,1],'or', linestyle='None',ms=1)
ax.set_aspect('equal')
plt.savefig('points_traslated.png',dpi=600)
plt.show()
#Perform relaxation
for i in range(self.relaxation_times):
field = Field(points)
field.relax()
points= field.get_points()
points=np.array(points)
fig, ax = plt.subplots(1, 1, figsize=(20,20))
ax.plot(points[:,0],points[:,1],'og', linestyle='None',ms=1)
ax.set_aspect('equal')
plt.savefig('points_traslated_relaxed.png',dpi=600)
plt.show()
points=self.de_translate(points)
points=np.array(points)
fig, ax = plt.subplots(1, 1, figsize=(20,20))
ax.plot(points[:,0],points[:,1],'om', linestyle='None',ms=1)
for k, dot in enumerate(self.polygon_external):
if k+1<len(self.polygon_external):
ax.plot([dot[0],self.polygon_external[k+1][0]],[dot[1],self.polygon_external[k+1][1]],'og')
ax.set_aspect('equal')
plt.savefig('points_relaxed.png',dpi=600)
plt.show()
#Drops points out of shapefile
rlx_points=[]
for point in points:
if self.check_in_out('external',[point[0],point[1]],'shapefile'):
rlx_points.append([point[0],point[1]])
rlx_points=np.array(rlx_points)
#Assign names to new points
position={}
for n, layer in enumerate(self.rock_dict):
for i,point in enumerate(rlx_points):
name=self.well_block_assign_from_relaxation(self.rock_dict[layer][0][0],i)
data['blocks'][name]=[n+1,point[0],point[1],self.rock_dict[layer][4],self.rock_dict[layer][5]]
if n==0:
position[i]=name
#Read border to clip write into in file
borders=gpd.read_file( self.outer_polygon)
border_points=[]
for line in borders.iterrows():
pointList = line[1].geometry.exterior.coords.xy
for point in zip(pointList[0],pointList[1]):
border_points.append([point[0],point[1]])
data['borders']=self.polygon_external
#Rewrite the wells dictionaries by finding the new closer position to the wells
wells_dictionary={}
well_corr={}
wells=pd.read_csv('../input/well_feedzone_xyz.csv',sep=',')
tree = cKDTree(rlx_points)
for index, well in wells.iterrows():
distance = tree.query([well['x'],well['y']])
block_assig=position[distance[1]]
well_corr[well['well']]=[block_assig[1:5]]
x= data['blocks'][block_assig][1]
y=data['blocks'][block_assig][2]
for n in range(self.number_of_layer):
block_assig_i=self.rock_dict[n+1][0]+block_assig[2:5]
wells_dictionary[block_assig_i]=[well['well'],\
x,\
y,\
well['type'],\
n+1,
data['blocks'][block_assig][3],\
self.rock_dict[n+1][2],\
data['blocks'][block_assig][4]]
json.dump(well_corr, open("../mesh/well_dict.txt",'w'),sort_keys=True, indent=4)
json.dump(wells_dictionary, open("../mesh/wells_correlative.txt",'w'),sort_keys=True, indent=4)
return data
[docs] def run_amesh_from_relaxation(self):
os.system("amesh\n")
files2move=['in','conne','eleme','segmt']
for fn in files2move:
shutil.move(fn,'../mesh/from_amesh/%s'%fn)
return None
[docs] def plot_voronoi(self):
"""En caso de ser solicitado, grafica la malla a partir de los archivo en la carpeta ../mesh/from_amesh
"""
welldict=self.well_blk_assign()['well_blocks_names']
if os.path.isfile('../mesh/from_amesh/segmt'):
data=np.genfromtxt('../mesh/from_amesh/segmt', dtype="f8,f8,f8,f8,i8,S5,S5", names=['X1','Y1','X2','Y2','index','elem1','elem2'],delimiter=[15,15,15,15,5,5,5])
fig, ax0 = plt.subplots(figsize=(10,10))
#Set of the principal plot
ax0.axis([self.Xmin,self.Xmax,self.Ymin,self.Ymax])
ax0.set_xlabel('East [m]')
ax0.set_ylabel('North [m]')
ax0.set_title("Mesh for layer %s"%self.layer_to_plot,y=1.04)
ax0.set_xticks(range(self.Xmin+self.x_from_boarder,self.Xmax+self.x_space-self.x_from_boarder,self.x_space))
ax0.set_yticks(range(self.Ymin+self.y_from_boarder,self.Ymax+self.y_space-self.y_from_boarder,self.y_space))
ax0.xaxis.set_minor_locator(AutoMinorLocator())
ax0.yaxis.set_minor_locator(AutoMinorLocator())
#Plot of the Yaxis in the top
ax1 = ax0.twinx()
ax1.set_ylim(ax0.get_ylim())
ax1.set_yticks(ax0.get_yticks())
ax1.yaxis.set_minor_locator(AutoMinorLocator())
#Plot of the Xaxis in the top
ax2 = ax0.twiny()
ax2.set_xticks(ax0.get_xticks())
ax2.set_xlim(ax0.get_xlim())
ax2.xaxis.set_minor_locator(AutoMinorLocator())
for n in np.unique(data['elem1']):
Xfill=[]
Yfill=[]
if n[0:2] in self.color_dict[self.layer_to_plot][0]:
for j in range(len(data['X1'])):
cnt=0
if data['elem1'][j]==n:
Xfill.append(data['X1'][j])
Xfill.append(data['X2'][j])
Yfill.append(data['Y1'][j])
Yfill.append(data['Y2'][j])
plt.plot([data['X1'][j],data['X2'][j]],[data['Y1'][j],data['Y2'][j]],'-k')
if self.plot_names:
if cnt==0:
if welldict[n][3]=='WELL':
ax0.text(welldict[n][1],welldict[n][2],\
"%s, %s"%(welldict[n][0],n),fontsize=8)
color_dot='r'
else:
ax0.text(welldict[n][1],welldict[n][2],\
welldict[n][0],fontsize=8)
color_dot='b'
cnt+=1
if self.plot_centers:
if welldict[n][3]=='WELL':
color_dot='r'
else:
color_dot='b'
ax0.plot(welldict[n][1],welldict[n][2],'%so'%color_dot,ms=1)
ax0.fill(Xfill,Yfill,welldict[n][6],alpha=0.5)
ax0.grid()
ax0.grid(which='major', color='#CCCCCC', linestyle='--',alpha=0.5)
fig.savefig('../mesh/plot_voronoi_'+str(datetime.datetime.now().strftime('%Y-%m-%d_%H%M%S')),bbox_inches="tight", pad_inches = 0)
plot_s=plt.show()
else:
plot_s="Excecute input_file_to_amesh() function first"
return plot_s
[docs] def to_steinar(self):
"""Convierte los archivos de salida amesh en formato de entrada para Steinar
"""
if os.path.isfile('../mesh/from_amesh/segmt'):
shutil.copyfile('../mesh/from_amesh/segmt', '../mesh/to_steinar/segmt')
ele_file_st=open('../mesh/to_steinar/eleme','w')
ele_file=open('../mesh/from_amesh/eleme','r')
data_eleme=ele_file.readlines()
if self.set_inac_from_poly or self.set_inac_from_inner:
source_file="../mesh/from_amesh/eleme"
data_eleme=pd.read_csv(source_file,delim_whitespace=True,skiprows=1,header=None,names=['block','rocktype','vol','X','Y','Z'])
data_eleme.set_index('block')
ele_file_st.write("eleme\n")
if self.set_inac_from_inner:
source='inner'
if self.set_inac_from_poly:
source='shapefile'
for index, row in data_eleme.iterrows():
check=self.check_in_out('internal',[row['X'],row['Y']],source=source)
outside=1
if not check:
outside=-1
ele_file_st.write("%5s%10s%5s%10.3E\n"%(row['block']," ",row['rocktype'],row['vol']*outside))
else:
cnt=0
ele_file_st.write("eleme")
for a_line in data_eleme:
if cnt!=0:
name_rock_vol = a_line[0:30]
ele_file_st.write("\n")
ele_file_st.write("%s"%name_rock_vol)
cnt+=1
ele_file_st.close()
ele_file.close()
conne_file_st=open('../mesh/to_steinar/conne','w')
conne_file=open('../mesh/from_amesh/conne','r')
data_conne=conne_file.readlines()
cnt=0
for a_line in data_conne:
if cnt!=0:
name_conne = a_line[0:11]
if '*' in name_conne:
pass
elif cnt<=len(data_conne)-2:
try:
conne_file_st.write("%s"%name_conne)
conne_file_st.write(" 0 0 0 ")
nod4=a_line[60:64]
if len(nod4)==4:
per=3
else:
per=a_line[29:30]
conne_file_st.write("%s"%per)
nod1=a_line[30:40].replace('+', '')
conne_file_st.write(" %9.3E"%(float(nod1)))
nod2=a_line[40:50].replace('+', '')
conne_file_st.write(" %9.3E"%(float(nod2)))
nod3=a_line[50:60].replace('+', '')
conne_file_st.write(" %9.3E"%(float(nod3)))
if len(nod4)==4:
#nod4=-1
conne_file_st.write("{:9.1f}".format(float(nod4)))
elif len(nod4)==1:
conne_file_st.write(" 0.0")
conne_file_st.write(" 0.0000E+00\n")
except ValueError:
pass
else:
conne_file_st.write("conne\n")
cnt+=1
conne_file_st.close()
conne_file.close()
in_file_st=open('../mesh/to_steinar/in','w')
in_file=open('../mesh/from_amesh/in','r')
data_in=in_file.readlines()
in_file_st.write("bound\n")
if not self.from_leapfrog:
Xarray=np.array([self.Xmin,self.Xmax])
Yarray=np.array([self.Ymin,self.Ymax])
else:
xmax,xmin,ymin,ymax,regular_mesh=self.from_leapfrog_mesh()
Xarray=np.array([xmin,xmax])
Yarray=np.array([ymin,ymax])
if self.rotate:
borders=self.polygon_external
for point in borders:
string_bound=" %6.2d %6.2d\n"%(point[0],point[1])
in_file_st.write(string_bound)
else:
in_file_st.write(" %6.2d %6.2d\n"%(Xarray[0],Yarray[0]))
in_file_st.write(" %6.2d %6.2d\n"%(Xarray[0],Yarray[1]))
in_file_st.write(" %6.2d %6.2d\n"%(Xarray[1],Yarray[1]))
in_file_st.write(" %6.2d %6.2d\n\n"%(Xarray[1],Yarray[0]))
in_file_st.write("\nlocat\n")
cnt=0
for ny in data_in:
if cnt>0 and cnt<(len(data_in)-10):
read1=ny[0:6]
read2=ny[8:11]
read3=ny[11:31]
read4=ny[31:52]
read5=ny[51:72]
read6=ny[71:92]
in_file_st.write("%s"%read1)
in_file_st.write("%4s"%read2)
in_file_st.write("%23.8E"%float(read3))
in_file_st.write("%24.8E"%float(read4))
in_file_st.write("%24.8E"%float(read5))
in_file_st.write("%24.8E\n"%float(read6))
cnt+=1
in_file_st.close()
in_file.close()
#Dummy rocktype
rock_file_st=open('../mesh/to_steinar/rocks','w')
rock_file_st.write("rock1 02.6500E+031.0000E-011.0000E-151.0000E-151.0000E-152.1000E+008.5000E+02 160 243 150")
rock_file_st.close()
else:
pass
[docs] def to_GIS(self):
"""Toma como entrada el archivo segment de ../mesh/from_amesh y eleme de ../mesh/to_stainar y los convierte en shapefile con atributos de roca, nombre y volumen
"""
if os.path.isfile('../mesh/from_amesh/segmt') and os.path.isfile('../mesh/to_steinar/eleme'):
if self.plot_all_GIS:
max_layer=self.number_of_layer
min_layer=1
else:
max_layer=self.layer_to_plot
min_layer=self.layer_to_plot
for ln in range(min_layer,max_layer+1,1):
w=shapefile.Writer('../mesh/GIS/mesh_layer_%s'%ln)
w.field('BLOCK_NAME', 'C', size=5)
w.field('ROCKTYPE', 'C', size=5)
w.field('VOLUMEN', 'F', decimal=10)
alttype = np.dtype([('X1', '<f8'), ('Y1', '<f8'), ('X2', '<f8'), ('Y2', '<f8'), ('index', '<f8'), ('elem1', 'U5'), ('elem2', 'U5')])
data = np.genfromtxt('../mesh/to_steinar/segmt', dtype=alttype, delimiter=[15,15,15,15,5,5,5])
elem_file = open('../mesh/to_steinar/eleme', 'r')
plain_elemfile=elem_file.read()
for n in np.unique(data['elem1']):
if n[0:2] in self.color_dict[ln][0]:
points=[]
line = re.findall(r"%s.*$"%n, plain_elemfile,re.MULTILINE)
rocktype=str(line)[17:22]
volumen=float(str(line)[22:32])
for j in range(len(data['X1'])):
if data['elem1'][j]==n:
point1=[data['X1'][j],data['Y1'][j]]
point2=[data['X2'][j],data['Y2'][j]]
points.append(point1)
points.append(point2)
w.poly([points])
w.record(n,rocktype,volumen)
w.close()
elem_file.close()
else:
None
[docs]def mesh_creation_func(input_mesh_dictionary,input_dictionary):
"""Creates a grid
Most of the parameters are keywords from the input dictionary input_mesh_dictionary. Otherwise, the input dictionary is specified.
Parameters
----------
filename : str
File name with well feedzone location
filepath : str
Path of input files
Xmin : float
Minimun X coordinates for the grid
Xmax : float
Maximun X coordinates for the grid
Ymin : float
Minimun Y coordinates for the grid
Ymax : float
Maximun Y coordinates for the grid
toler : float
AMESH parameter
layers : dictionary
Name (correlative) and thickness of every layer on the model, keyword on input_dictionary
layer_to_plot : int
In case it is specified a voronoi plot will be performed
x_space : float
Horizontal distance between elements for the outerfield
y_space : float
Vertical distance between elements for the outerfield
radius_criteria: float
Minimun distance between well location and a regular element
x_from_boarder: float
Horizontal distance from the first element to the east border
y_from_boarder: float
Vertical distance from the first element to the south border
x_gap_min: float
Minimun X coordinates on the grid for the well field
x_gap_max: float
Maximun X coordinates on the grid for the well field
x_gap_space: float
Horizontal distance between elements for the farfield
y_gap_min: float
Minimun Y coordinates on the grid for the well field o
y_gap_max: float
Maximun X coordinates on the grid for the well field
y_gap_space: float
Vertical distance between elements for the farfield
plot_names: bool
If true it plots the name of the blocks from the selected layer to plot
plot_centers: bool
If true it plots the centers of the blocks from the selected layer to plot
z0_level: float
Reference level (elevation) for all the grid, keyword on input_dictionary
mesh_creation: bool
If true the mesh is created
plot_layer: bool
If true it plots the selected layer
to_steinar: bool
If true it creates the input files for steinar
to_GIS: bool
If true it generates a shapefile of the selected layer
plot_all_GIS: bool
If true it generates a shapefile of all layers
from_leapfrog: bool
lee archivos leapfrong ../mesh/from_leapfrog/LF_geometry.dat y ../mesh/from_leapfrog/LF_t2.dat, sin embargo se pierde la simbologia usada en leapfrog y no te utiliza la malla regular ni los pozos. Solamente se crea la malla utilizando amesh
line_file: str
It defines the path and name of a line that can represented a fault or other structure on the mesh, The input file must contain the header: ID,X,Y on csv format. ID referes to the same
structure, thus, more than one structure can be defined on a single file.
fault_distance: float
In case a line_file is define, some paralels elements will be created at a defined distance
with_polygon: bool
If true a shapefile will be read to define the wellfield.
polygon_shape: str
The shapefile deines the wellfield boundaries. The shape must not contain any cavity
set_inac_from_poly: bool
If true all the elements on the outside of the shapefile are defined as inactive
set_inac_from_inner:bool
If true all the elements on the outerfield are defined as inactive
rotate: bool
If true it rotates the mesh a defined angle
angle: float
Angle in degrees
inner_mesh_type: string
Type of mesh on the inner part of the mesh, it could be 'honeycomb' or 'regular'
Returns
-------
file
eleme: list of blocks from the grid
file
conne : list of connections on the grid
shapefile
mesh_{field}_layer_{layer} : shapefile of a defined (or all) layer including rock distribution
plot
Voronoi plot (in case is specified)
Attention
---------
A copy of AMESH must be on the path or directory
"""
layers_thick=list(map(float,np.array(list(input_dictionary['LAYERS'].values()))[:][:,1]))
blocks=py2amesh(input_mesh_dictionary['filename'],input_mesh_dictionary['filepath'],input_mesh_dictionary['Xmin'],input_mesh_dictionary['Xmax'],input_mesh_dictionary['Ymin'],input_mesh_dictionary['Ymax'],\
input_mesh_dictionary['toler'],layers_thick,input_mesh_dictionary['layer_to_plot'],input_mesh_dictionary['x_space'],input_mesh_dictionary['y_space'],input_mesh_dictionary['radius_criteria'],input_mesh_dictionary['x_from_boarder'],input_mesh_dictionary['y_from_boarder'],\
input_mesh_dictionary['x_gap_min'],input_mesh_dictionary['x_gap_max'],input_mesh_dictionary['x_gap_space'],input_mesh_dictionary['y_gap_min'],input_mesh_dictionary['y_gap_max'],input_mesh_dictionary['y_gap_space'],input_mesh_dictionary['plot_names'],input_mesh_dictionary['plot_centers'],\
input_dictionary['z_ref'],input_mesh_dictionary['plot_all_GIS'],input_mesh_dictionary['from_leapfrog'],input_mesh_dictionary['line_file'],input_mesh_dictionary['fault_distance'],input_mesh_dictionary['with_polygon'],input_mesh_dictionary['polygon_shape'],\
input_mesh_dictionary['set_inac_from_poly'],input_mesh_dictionary['set_inac_from_inner'],input_mesh_dictionary['rotate'],input_mesh_dictionary['angle'],input_mesh_dictionary['inner_mesh_type'],\
input_mesh_dictionary['distance_points'],input_mesh_dictionary['fault_rows'],input_mesh_dictionary['relaxation_times'],input_mesh_dictionary['points_around_well'],input_mesh_dictionary['distance_points_around_well'], input_mesh_dictionary['outer_polygon'])
#data=blocks.relaxation()
if input_mesh_dictionary['mesh_creation']:
if input_mesh_dictionary['relaxation_times']>0:
blocks.input_file_to_amesh_from_relaxation()
blocks.run_amesh_from_relaxation()
else:
blocks.input_file_to_amesh()
if input_mesh_dictionary['to_steinar']:
blocks.to_steinar()
if input_mesh_dictionary['plot_layer']:
blocks.plot_voronoi()
if input_mesh_dictionary['to_GIS']:
blocks.to_GIS()
return None #data #None
[docs]def change_ref_elevation(variation=0):
"""It modifies the in files from the folders to_steinar and from_amesh by increasing (or decreasing) a fixed value
Parameters
----------
variation: float
Defined value change the elevation of the in file from amesh
Returns
-------
file
in: modified on folders ../mesh/to_steinar and ../mesh/from_amesh
Attention
---------
The change is generated on all the elements. It is recommened to run ELEM_to_json() after this execution
Examples
--------
>>> change_ref_elevation(variation=-100)
"""
in_file_steinar='../mesh/to_steinar/in'
in_file_amesh='../mesh/from_amesh/in'
if os.path.isfile(in_file_steinar) and os.path.isfile(in_file_amesh):
steinar_colums=[(0,5),(5,10),(10,35),(35,60),(60,85),(85,110)]
steinar_data=pd.read_fwf(in_file_steinar,colspecs=steinar_colums,skiprows=7,header=None,names=['block','level','X','Y','Z','h'])
steinar_data['Z']=steinar_data['Z'] + variation
string=""
with open(in_file_steinar,'r') as f:
for line in f.readlines():
string+=line
if 'locat' in line:
break
for index, row in steinar_data.iterrows():
string+="{:5s}{:5d}{:23.8f}{:24.8f}{:24.8f}{:24.8f}\n".format(row['block'], row['level'],row['X'],row['Y'],row['Z'],row['h'])
steinar_in_file=open(in_file_steinar,'w')
steinar_in_file.write(string)
steinar_in_file.close()
amesh_colums=[(0,5),(5,10),(10,30),(30,50),(50,70),(70,90)]
amesh_data=pd.read_fwf(in_file_amesh,colspecs=amesh_colums,skiprows=1,header=None,names=['block','level','X','Y','Z','h'])
amesh_data.dropna(inplace=True)
amesh_data['Z']=amesh_data['Z'] + variation
string="locat\n"
for index, row in amesh_data.iterrows():
string+="{:5s}{:5.0f}{:20.3f}{:20.3f}{:20.3f}{:20.3f}\n".format(row['block'], row['level'],row['X'],row['Y'],row['Z'],row['h'])
with open(in_file_amesh,'r') as f:
check_in=False
for line in f.readlines():
if 'bound' in line:
check_in=True
string+='\n'
if check_in:
string+=line
amesh_in_file=open(in_file_amesh,'w')
amesh_in_file.write(string)
amesh_in_file.close()
[docs]def empty_mesh():
"""If the suggested T2GEORES folder structure is used. The folder mesh/to_steinar, mesh/from_amesh and mesh/GIS are emptied
Note
----
No input parameter is needed
Examples
--------
>>> empty_mesh()
"""
folders=['../mesh/to_steinar','../mesh/from_amesh','../mesh/GIS']
for folder in folders:
for file in os.listdir(folder):
try:
file_path=os.path.join(folder,file)
os.remove(file_path)
except IsADirectoryError:
pass
[docs]def ELEM_to_json(input_dictionary,to_sql=False):
"""It combines the files eleme and in from the steinar folder into a single json file
Parameters
----------
input_dictionary: dictionary
Dictionary containing the path and name of database and the path of the input file
Returns
-------
file
ELEME.json: on ../mesh/
Attention
---------
Currently, it uses the steinar function because it assume the modeler will perfome the block assignation on it
Examples
--------
>>> ELEM_to_json()
"""
source_file="../mesh/to_steinar/in"
eleme_file="../mesh/to_steinar/eleme"
eleme_dict={}
col_eleme=[(0,6),(15,20),(20,30)]
if os.path.isfile(source_file) and os.path.isfile(eleme_file):
data_in=pd.read_csv(source_file,delim_whitespace=True,skiprows=7,header=None,names=['ELEME','LAYER_N','X','Y','Z','h'])
data_eleme=pd.read_fwf(eleme_file,colspecs=col_eleme,skiprows=1,header=None,names=['ELEME','MA1','VOLX'])
#data_eleme[['rocktype','vol']] =data_eleme.data.apply(lambda x: pd.Series([str(x)[0:5],str(x)[5:16]]))
#data_eleme.drop('data', axis=1, inplace=True)
data_in.set_index('ELEME')
data_eleme.set_index('ELEME')
data_eleme=data_eleme.merge(data_in, left_on='ELEME', right_on='ELEME',how='left')
data_eleme.set_index('ELEME',inplace=True)
if to_sql:
data_eleme['model_version']=[input_dictionary['model_version']]*len(data_eleme)
data_eleme['model_timestamp']=[datetime.datetime.now()]*len(data_eleme)
conn=sqlite3.connect(input_dictionary['db_path'])
data_eleme.to_sql('ELEME',if_exists='append',con=conn,index=True)
conn.close()
data_eleme.to_json('../mesh/ELEME.json',orient="index",indent=2)
else:
sys.exit("The file %s or directory do not exist"%source_file)
[docs]def CONNE_to_json(input_dictionary,to_sql=False):
"""It creates a json file from the CONNE file on mesh/to_steinar folder
Parameters
----------
input_dictionary: dictionary
Dictionary containing the path and name of database and the path of the input file
Returns
-------
file
CONNE.json: on ../mesh/
Attention
---------
Currently, it uses the steinar folder because it assume the modeler will perfome the block assignation on it
Examples
--------
>>> CONNE_to_json()
"""
source_file="../mesh/to_steinar/conne"
conne_dict={}
if os.path.isfile(source_file):
with open(source_file) as rock_file:
for line in rock_file.readlines()[1:]:
linex=line.split()
if float(linex[8])<=0:
conne_dict[linex[0]]=[linex[0][0:5],linex[0][5:10],int(linex[4]),float(linex[5]),float(linex[6]),float(linex[7]),float(linex[8])]
else:
conne_dict[linex[0]]=[linex[0][0:5],linex[0][5:10],int(linex[4]),float(linex[5]),float(linex[6]),float(linex[7]),-1*float(linex[8])]
CONNE_pd=pd.DataFrame.from_dict(conne_dict,orient='index', columns=['ELEME1','ELEME2','ISOT','D1','D2','AREAX', 'BETAX'])
CONNE_pd.to_json('../mesh/CONNE.json',orient="index",indent=2)
if to_sql:
CONNE_pd['model_version']=[input_dictionary['model_version']]*len(CONNE_pd)
CONNE_pd['model_timestamp']=[datetime.datetime.now()]*len(CONNE_pd)
conn=sqlite3.connect(input_dictionary['db_path'])
CONNE_pd.to_sql('CONNE',if_exists='append',con=conn,index=False)
conn.close()
else:
sys.exit("The file %s or directory do not exist"%source_file)
[docs]def segmnt_to_json(input_dictionary,to_sql=False):
"""It combines the information from eleme and segmt from folder mesh/to_steinar into a single json
Parameters
----------
input_dictionary: dictionary
Dictionary containing the path and name of database and the path of the input file
Returns
-------
file
segmt.json: on ../mesh/
Attention
---------
Currently, it uses the steinar folder because it assume the modeler will perfome the block assignation on it.
Examples
--------
>>> segmnt_to_json()
"""
eleme_file="../mesh/to_steinar/eleme"
segmt_file="../mesh/to_steinar/segmt"
elment_dict={}
col_eleme=[(0,6),(15,20),(20,30)]
col_segment=[(0,15),(15,30),(30,45),(45,60),(60,65),(65,70),(70,75)]
if os.path.isfile(segmt_file) and os.path.isfile(eleme_file):
segmt_data=pd.read_fwf(segmt_file,colspecs=col_segment,header=None,names=['X1','Y1','X2','Y2','redundant','ELEME1','ELEME2'])
if to_sql:
segmt_data['model_version']=[input_dictionary['model_version']]*len(segmt_data)
segmt_data['model_timestamp']=[datetime.datetime.now()]*len(segmt_data)
conn=sqlite3.connect(input_dictionary['db_path'])
segmt_data.to_sql('segment',if_exists='append',con=conn,index=False)
conn.close()
eleme_data=pd.read_fwf(eleme_file,colspecs=col_eleme,skiprows=1,header=None,names=['block','rocktype','vol'])
eleme_data=eleme_data.merge(segmt_data, left_on='block', right_on='ELEME1',how='left')
eleme_data['points'] = eleme_data[['X1','Y1','X2','Y2']].apply(lambda r: tuple(r), axis=1).apply(np.array)
eleme_data=eleme_data.groupby(['block','rocktype','vol'])['points'].agg(lambda x: list(np.concatenate(x.values))).reset_index()
eleme_data.set_index('block', inplace=True)
eleme_data.to_json('../mesh/segmt.json',orient="index",indent=2)
else:
sys.exit("Either the file %s or %s do not exist"%(eleme_file,segmt_file))
[docs]def plot_voronoi(input_mesh_dictionary, geners,layer='D',plot_center=False,mark_elements=False,cross_section=None,savefig=False):
"""It plots a selected layer, showing the rock distribution
Returns
-------
layer : str
It defines the layer correlative to plot
plot_center : bool
If True it plots the center point for each element
mark_elements : bool
If True it considers the final four identifiers from each constant sink/source element defined on geners and fill the block with a different color.
cross_section : bool
If true plot a line where a vertical cross section is created with the output module
savefig : bool
If true a png file is saved on ../mesh
input_mesh_dictionary : dictionary
Contains a unique color for every rocktype.
geners : dictionary
Contains the neccesary information to define a sink/source. g.e. 'DA110':{'SL':'GEN','NS':10,'TYPE':'MASS','GX':1,'EX':1.1E6}
Returns
-------
plot
layer_{layer}: containing the rock distribution from the selected layer
Attention
---------
The functions ELEM_to_json()
Examples
--------
>>> segmnt_to_json() and segmnt_to_json() should be executed previously
"""
segmt_json_file='../mesh/segmt.json'
eleme_json_file='../mesh/ELEME.json'
wells_corr_json_file='../mesh/wells_correlative.txt'
if cross_section:
x_points=cross_section['cross_section']['x']
y_points=cross_section['cross_section']['y']
if os.path.isfile(segmt_json_file):
with open(segmt_json_file) as file:
blocks=json.load(file)
rocktypes={}
cnt=1
for block in blocks:
if blocks[block]['rocktype'] not in rocktypes:
rocktypes[blocks[block]['rocktype']]=input_mesh_dictionary['colors'][cnt]
cnt+=1
if os.path.isfile(eleme_json_file):
with open(eleme_json_file) as file:
ELEME=json.load(file)
if mark_elements:
legend='Legend'
if os.path.isfile(wells_corr_json_file):
with open(wells_corr_json_file) as file:
well_corr=json.load(file)
for gen in geners:
src=geners[gen]['SL']+str(geners[gen]['NS'])
corr=gen[1:6]
well_corr[src]=[corr]
corr_well = {y[0]:x for x,y in well_corr.items()}
else:
legend='ROCKS'
fig, ax0 = plt.subplots(figsize=(10,10))
ax0.set_xlabel("East [m]")
ax0.set_ylabel("North [m]")
if cross_section:
ax0.plot(x_points,y_points,linestyle='--',color='darkred')
ax0.text(x_points[0], y_points[0], 'A',color='darkred',horizontalalignment='left')
ax0.text(x_points[-1], y_points[-1], "A'",color='darkred',horizontalalignment='left')
colors_on_plot=[]
legends=[]
for block in blocks:
if block[0]==layer:
for i, point in enumerate(blocks[block]['points']):
if i==0:
fill_x=[]
fill_y=[]
line=[]
if i!=0 and (i)%4==0:
plt.plot(line[0],line[1],line[2],line[3],'-k')
fill_x.extend([line[0],line[2]])
fill_y.extend([line[1],line[3]])
line=[]
if i%4 in [0,1,2,3]:
line.append(point)
if plot_center:
ax0.plot(ELEME[block]['X'],ELEME[block]['Y'],'.k')
ax0.fill(fill_x,fill_y,facecolor=rocktypes[blocks[block]['rocktype']], edgecolor='black',alpha=0.5)
if mark_elements:
if 'red' not in colors_on_plot:
legend_for_rock = mpatches.Patch(color='red', label='sinks/sources',alpha=0.75)
colors_on_plot.append('red')
legends.append(legend_for_rock)
if any(corr[0] in block for corr in well_corr.values()):
ax0.fill(fill_x,fill_y,facecolor='red', edgecolor='black',alpha=0.75)
ax0.annotate(corr_well[block[1:6]],
xy=(ELEME[block]['X'],ELEME[block]['Y']), xycoords='data',
xytext=(ELEME[block]['X']+500,ELEME[block]['Y']+500),textcoords='data',fontsize=10,
arrowprops=dict(arrowstyle="->", color="white",
shrinkA=5, shrinkB=5,
patchA=None, patchB=None,
connectionstyle="bar,angle=180,fraction=-0.2"))
if rocktypes[blocks[block]['rocktype']] not in colors_on_plot:
legend_for_rock = mpatches.Patch(color=rocktypes[blocks[block]['rocktype']], label=blocks[block]['rocktype'],alpha=0.5)
colors_on_plot.append(rocktypes[blocks[block]['rocktype']])
legends.append(legend_for_rock)
plt.subplots_adjust(right=0.8)
plt.legend(handles=legends,title=legend, bbox_to_anchor=(1.01, 1), loc='upper left')
plt.axis('equal')
if savefig:
plt.savefig('../mesh/layer_%s.png'%layer,dpi=300,transparent=False,bbox_inches='tight')
plt.show()
else:
sys.exit("The file %s does not exist"%(segmt_json_file))
[docs]def to_GIS(input_mesh_dictionary, layer='A'):
"""It plots a selected layer, showing the rock distribution
Returns
-------
layer : str
It defines the layer correlative to plot
plot_center : bool
If True it plots the center point for each element
mark_elements : bool
If True it considers the final four identifiers from each constant sink/source element defined on geners and fill the block with a different color.
cross_section : bool
If true plot a line where a vertical cross section is created with the output module
savefig : bool
If true a png file is saved on ../mesh
input_mesh_dictionary : dictionary
Contains a unique color for every rocktype.
geners : dictionary
Contains the neccesary information to define a sink/source. g.e. 'DA110':{'SL':'GEN','NS':10,'TYPE':'MASS','GX':1,'EX':1.1E6}
Returns
-------
plot
layer_{layer}: containing the rock distribution from the selected layer
Attention
---------
The functions ELEM_to_json()
Examples
--------
>>> segmnt_to_json() and segmnt_to_json() should be executed previously
"""
segmt_json_file='../mesh/segmt.json'
eleme_json_file='../mesh/ELEME.json'
wells_corr_json_file='../mesh/wells_correlative.txt'
if os.path.isfile(segmt_json_file):
with open(segmt_json_file) as file:
blocks=json.load(file)
if os.path.isfile(eleme_json_file):
with open(eleme_json_file) as file:
ELEME=json.load(file)
w=shapefile.Writer('../mesh/GIS/mesh_layer_%s'%layer)
w.field('BLOCK_NAME', 'C', size=5)
w.field('ROCKTYPE', 'C', size=5)
w.field('VOLUMEN', 'F', decimal=10)
colors_on_plot=[]
legends=[]
for block in blocks:
if block[0]==layer:
points=[]
for i, point in enumerate(blocks[block]['points']):
if i%2==0:
x=point
else:
y=point
points.append([x,y])
w.poly([points])
w.record(block,blocks[block]['rocktype'],blocks[block]['vol'])
w.close()
else:
sys.exit("The file %s does not exist"%(segmt_json_file))
[docs]def CONNE_count(save_csv=False):
"""
Return a dataframe with the name of elements and the number of connection
Parameters
----------
save_csv: bool
If true writes the output dataframe as csv, default False
Returns
-------
dataframe
count_conne: blocks and number of connections
file:
CONNE_count.csv: at ../mesh/
"""
conne_file='../mesh/CONNE.json'
if os.path.isfile(conne_file):
CONNE_df=pd.read_json(conne_file)
CONNE_df=CONNE_df.T
CONNE_df['BLOCKS']=CONNE_df.index
CONNE_df[['block1','block2']] =CONNE_df.BLOCKS.apply(lambda x: pd.Series([str(x)[0:5],str(x)[5:10]]))
count_conne_1=pd.DataFrame(CONNE_df['block1'].value_counts().reset_index())
count_conne_2=pd.DataFrame(CONNE_df['block2'].value_counts().reset_index())
count_conne=pd.concat([count_conne_1,count_conne_2])
count_conne=count_conne.groupby('index').sum()
count_conne['total']=count_conne['block1']+count_conne['block2']
del count_conne['block1']
del count_conne['block2']
count_conne=count_conne.sort_values(by=['total'],ascending=False)
if save_csv:
count_conne.to_csv('../mesh/CONNE.csv')
return count_conne
else:
sys.exit("The file %s or directory do not exist"%conne_file)
[docs]def mesh_to_paraview(input_dictionary):
"""
It generates a vtu file to be read on Paraview including the rocktype and block name.
Parameters
----------
input_dictionary: dictionary
Contains the information of the layers on the model.
Returns
-------
file
model_mesh.vtu: on ../mesh
"""
segmt_json_file='../mesh/segmt.json'
if os.path.isfile(segmt_json_file):
with open(segmt_json_file) as file:
blocks=json.load(file)
else:
sys.exit("The file %s or directory do not exist, run segmnt_to_json from regeo_mesh"%segmt_json_file)
ugrid=vtk.vtkUnstructuredGrid()
block_name=vtk.vtkStringArray()
block_name.SetName('block')
rocktype=vtk.vtkStringArray()
rocktype.SetName('rocktype')
layers_dict=geometry.vertical_layers(input_dictionary)
layers_name={}
for j, layer in enumerate(layers_dict['name']):
layers_name[layer]=[layers_dict['top'][j],layers_dict['bottom'][j]]
points_vtk=vtk.vtkPoints()
#elem_file='../mesh/vtu/ELEME.json'
elem_file='../mesh/ELEME.json'
if os.path.isfile(elem_file):
with open(elem_file) as file:
blocks_eleme=json.load(file)
else:
sys.exit("The file %s or directory do not exist, run segmnt_to_json from regeo_mesh"%elem_file)
cnt=0
for block in blocks:
if block in blocks_eleme.keys():
faceId=vtk.vtkIdList()
block_name.InsertNextValue(block)
rocktype.InsertNextValue(blocks[block]['rocktype'])
points=[]
for i, point in enumerate(blocks[block]['points']):
if i%2==0:
point_x=point
elif i%2==1:
point_y=point
if i%2==1 and i!=0 and [point_x,point_y] not in points:
points.append([point_x,point_y])
points=np.array(points)
hull = ConvexHull(points)
dict_points={}
cnt_int=0
for n in range(len(points[hull.vertices,0])):
dict_points[n]=[points[hull.vertices,0][::-1][n],points[hull.vertices,1][::-1][n],layers_name[block[0]][0]]
cnt_int+=1
for n in range(len(points[hull.vertices,0])):
dict_points[n+len(points)]=[points[hull.vertices,0][::-1][n],points[hull.vertices,1][::-1][n],layers_name[block[0]][1]]
cnt_int+=1
for key in dict_points:
points_vtk.InsertNextPoint(dict_points[key])
faceId.InsertNextId(len(points)+2)
faces=[[] for n in range(len(points)+2)]
faces[0]=[face+cnt for face in range(len(points))]
faces[1]=[face+len(points)+cnt for face in range(len(points))]
for nx in range(len(points)):
if nx+1<len(points):
faces[nx+2]=[nx+cnt,nx+len(points)+cnt,nx+1+len(points)+cnt,nx+1+cnt]
else:
faces[nx+2]=[nx+cnt,nx+len(points)+cnt,nx+cnt+1,cnt]
cnt+=cnt_int
for face in faces:
faceId.InsertNextId(len(face))
[faceId.InsertNextId(i) for i in face]
ugrid.InsertNextCell(vtk.VTK_POLYHEDRON,faceId)
ugrid.GetCellData().AddArray(block_name)
ugrid.GetCellData().AddArray(rocktype)
ugrid.SetPoints(points_vtk)
writer=vtk.vtkXMLUnstructuredGridWriter()
writer.SetInputData(ugrid)
writer.SetFileName('../mesh/vtu/model_mesh.vtu')
writer.SetDataModeToAscii()
writer.Update()
writer.Write()
[docs]def surface_cut(input_mesh_dictionary):
"""It generates an json file with the elements below the surface
Parameters
----------
input_mesh_dictionary: dictionary
Contains the path for the file on csv format with the topography data
Returns
-------
file
ELEME.json: on the path ../mesh/vtu/
Attention
---------
csv file should have X,Y,Z on the first line and comma (,) as a delimiter
"""
data_file=input_mesh_dictionary['surface_data']
surface=pd.read_csv(data_file,index_col=None,delimiter=';')
X=list(surface['X'])
Y=list(surface['Y'])
Z=list(surface['Z'])
elem_file='../mesh/ELEME.json'
if os.path.isfile(elem_file):
with open(elem_file) as file:
blocks=json.load(file)
X_eleme=[]
Y_eleme=[]
for block in blocks:
if block[0]=='A':
X_eleme.append(blocks[block]['X'])
Y_eleme.append(blocks[block]['Y'])
Zi=griddata((X,Y),Z,(X_eleme,Y_eleme),method='linear')
eleme_suf=np.array((X_eleme,Y_eleme,Zi)).T
tree=cKDTree(eleme_suf)
eleme_suf_limit={}
for i,block in enumerate(blocks):
if block[0]=='A':
eleme_suf_limit[block[1:6]]=Zi[i]
else:
break
delete=[]
for i,block in enumerate(blocks):
blocks[block]['OUT']=False
eq_suf_elevation_i=tree.query([blocks[block]['X'],blocks[block]['Y'],blocks[block]['Z']])[1]
if blocks[block]['Z']>eleme_suf[eq_suf_elevation_i][2]:
blocks[block]['OUT']=True
delete.append(block)
copy_eleme=blocks
for key in delete: del copy_eleme[key]
json.dump(copy_eleme, open("../mesh/vtu/ELEME.json",'w'),sort_keys=True, indent=4)
[docs]def clip_steinar_data():
"""Modifies the eleme file on steinar folder. It converts all elements on top of topography to ATMOS rock type and add rocktype to rock file
Attention
---------
The output file is written on a file named rocks2
"""
copy_eleme_file="../mesh/vtu/ELEME.json"
with open(copy_eleme_file) as file:
copy_eleme=json.load(file)
eleme_file_st=open('../mesh/to_steinar/eleme2','w')
eleme_file=open('../mesh/to_steinar/eleme','r')
data_eleme=eleme_file.readlines()
cnt=0
for line in data_eleme:
if (line[0:5]) in copy_eleme.keys():
eleme_file_st.write(line)
elif line[0:5]!='eleme':
eleme_file_st.write(line[0:5]+" "+'ATMOS'+line[20:-1]+'\n')
eleme_file_st.close()
eleme_file.close()
rock_file=open('../mesh/to_steinar/rocks','r')
rock_file2=open('../mesh/to_steinar/rocks2','w')
rock_data=rock_file.readlines()
cnt=0
for line in rock_data:
rock_file2.write(line)
rock_file2.write("ATMOS 02.6500E+030.9990E+001.0000E-121.0000E-121.0000E-122.1000E+008.5000E+02 0 0 0")
rock_file2.close()
rock_file.close()
[docs]def reasig_feedzone(input_dictionary):
"""It reassing the block related with a well feedzone
Attention
---------
Files: '../input/well_feedzone_xyz.csv' and '../mesh/ELEME.json' need to be updated
Returns
-------
file:
well_dict.txt: at ../mesh/
file:
wells_correlative.txt: at ../mesh/
"""
layers_info=geometry.vertical_layers(input_dictionary)
layer_middle={layers_info['name'][n]:layers_info['middle'][n] for n in range(len(layers_info['name']))}
layers_thickness = {input_dictionary['LAYERS'][k][0] : input_dictionary['LAYERS'][k][1] for k in input_dictionary['LAYERS']}
elem_file='../mesh/ELEME.json'
if os.path.isfile(elem_file):
with open(elem_file) as file:
blocks=json.load(file)
rlx_points = []
block_name = []
for block in blocks:
if block[0]=='A':
rlx_points.append([blocks[block]['X'],blocks[block]['Y']])
block_name.append(block)
wells_dictionary={}
well_corr={}
wells=pd.read_csv('../input/well_feedzone_xyz.csv',sep=',')
tree = cKDTree(rlx_points)
for index, well in wells.iterrows():
distance = tree.query([well['x'],well['y']])
print(distance)
block_assig=block_name[distance[1]]
well_corr[well['well']]=[block_assig[1:5]]
x= blocks[block_assig]['X']
y=blocks[block_assig]['Y']
for n, layer in enumerate(layer_middle):
block_assig_i=layer+block_assig[1:5]
wells_dictionary[block_assig_i]=[well['well'],\
x,\
y,\
well['type'],\
n+1,
layer_middle[block_assig[0]],\
"red",\
layers_thickness[block_assig[0]]]
json.dump(well_corr, open("../mesh/well_dict.txt",'w'),sort_keys=True, indent=4)
json.dump(wells_dictionary, open("../mesh/wells_correlative.txt",'w'),sort_keys=True, indent=4)
[docs]def set_layer_inactive(layers = []):
"""It sets innactive one layer from the mesh on the ELEME file
Parameters
----------
layer: str
Layer to set inactive
Returns
-------
file
in: modified on folders ../mesh/to_steinar and ../mesh/from_amesh
Examples
--------
>>> set_layer_inactive(variation = 'A')
"""
source_file="../mesh/to_steinar/in"
eleme_file="../mesh/to_steinar/eleme"
eleme_dict={}
col_eleme=[(0,6),(15,20),(20,30)]
if os.path.isfile(source_file) and os.path.isfile(eleme_file):
data_in=pd.read_csv(source_file,delim_whitespace=True,skiprows=7,header=None,names=['ELEME','LAYER_N','X','Y','Z','h'])
data_eleme=pd.read_fwf(eleme_file,colspecs=col_eleme,skiprows=1,header=None,names=['ELEME','MA1','VOLX'])
for layer in layers:
data_eleme.loc[(data_eleme['ELEME'].str[0] == layer) & (data_eleme.VOLX > 0), "VOLX"] = -data_eleme.VOLX
data_in.set_index('ELEME')
data_eleme.set_index('ELEME')
data_eleme=data_eleme.merge(data_in, left_on='ELEME', right_on='ELEME',how='left')
data_eleme.set_index('ELEME',inplace=True)
data_eleme.to_json('../mesh/ELEME.json',orient="index",indent=2)
else:
sys.exit("The file %s or directory do not exist"%source_file)
[docs]def tracks_to_paraview(input_dictionary):
"""It generates a line for each well track as vtu
Parameters
----------
input_dictionary: dictionary
List of well to be included
Returns
-------
file
tracks_{well}.vtu in input/survey/vtu/
"""
from vtkmodules.vtkCommonDataModel import (
vtkCellArray,
vtkPolyData,
vtkPolyLine
)
from vtkmodules.vtkCommonCore import vtkPoints
wells=[]
for key in ['WELLS','MAKE_UP_WELLS','NOT_PRODUCING_WELL']:
try:
for well in input_dictionary[key]:
wells.append(well)
except KeyError:
pass
conn=sqlite3.connect(input_dictionary['db_path'])
c=conn.cursor()
for well in wells:
data=pd.read_sql_query("SELECT MeasuredDepth FROM survey WHERE well='%s' ORDER BY MeasuredDepth;"%well,conn)
x_real=[]
y_real=[]
z_real=[]
well_proj=[]
data_xyz = {'x':[],
'y':[],
'z':[]}
for v in range(len(data['MeasuredDepth'])):
xf,yf,zf=geometry.MD_to_TVD(well,data['MeasuredDepth'][v])
data_xyz['x'].append(float(xf))
data_xyz['y'].append(float(yf))
data_xyz['z'].append(float(zf))
# Create a vtkPoints object and store the points in it
points = vtkPoints()
# Create a cell array to store the lines in and add the lines to it
cells = vtkCellArray()
# Create a polydata to store everything in
polyData = vtkPolyData()
for i, k in enumerate(data_xyz['x']):
p = [data_xyz['x'][i],data_xyz['y'][i],data_xyz['z'][i]]
points.InsertNextPoint(p)
n_point = len(data_xyz['x'])
polyLine = vtkPolyLine()
polyLine.GetPointIds().SetNumberOfIds(n_point)
for i in range(0, n_point):
polyLine.GetPointIds().SetId(i, i)
cells.InsertNextCell(polyLine)
# Add the points to the dataset
polyData.SetPoints(points)
# Add the lines to the dataset
polyData.SetLines(cells)
ugrid=vtk.vtkUnstructuredGrid()
ugrid.SetCells(vtk.VTK_POLY_LINE,cells)
ugrid.SetPoints(points)
writer=vtk.vtkXMLUnstructuredGridWriter()
writer.SetInputData(ugrid)
writer.SetFileName('../input/survey/vtu/tracks_%s.vtu'%well)
writer.SetDataModeToAscii()
writer.Update()
writer.Write()
[docs]def feedpoints_to_vtu():
"""It generates a cube (Hexahedron) for each feedzone
Returns
-------
file
{index}_{well}.vtu in output/vtu/feedzones/
Attention
---------
The file input/well_feedzone_xyz.csv has to be created prior the execution of this function
"""
from vtkmodules.vtkCommonDataModel import (
vtkCellArray,
vtkHexahedron
)
from vtkmodules.vtkCommonCore import vtkPoints
feedzones=pd.read_csv('../input/well_feedzone_xyz.csv',delimiter=",")
for index,row in feedzones.iterrows():
x = row['x']
y = row['y']
z = row['z']
well = row['well']
r = 25
xs = r*2**0.5
zs = r
points = vtkPoints()
numberOfVertices = 8
points = vtkPoints()
block_name=vtk.vtkStringArray()
block_name.SetName('block')
block_name.InsertNextValue('%s_%s'%(well,index))
points.InsertNextPoint(x-xs, y-xs, z-zs) # Face 1
points.InsertNextPoint(x+xs, y-xs, z-zs)
points.InsertNextPoint(x+xs, y+xs, z-zs)
points.InsertNextPoint(x-xs, y+xs, z-zs)
points.InsertNextPoint(x-xs, y-xs, z+zs) # Face 1
points.InsertNextPoint(x+xs, y-xs, z+zs)
points.InsertNextPoint(x+xs, y+xs, z+zs)
points.InsertNextPoint(x-xs, y+xs, z+zs)
# Create a hexahedron from the points
hexhedr = vtkHexahedron()
for i in range(0, numberOfVertices):
hexhedr.GetPointIds().SetId(i, i)
# Add the points and hexahedron to an unstructured grid
uGrid = vtk.vtkUnstructuredGrid()
uGrid.SetPoints(points)
uGrid.InsertNextCell(hexhedr.GetCellType(), hexhedr.GetPointIds())
uGrid.GetCellData().AddArray(block_name)
writer=vtk.vtkXMLUnstructuredGridWriter()
writer.SetInputData(uGrid)
writer.SetFileName('../output/vtu/feedzones/%s_%s.vtu'%(well,index))
writer.SetDataModeToAscii()
writer.Update()
writer.Write()