Import heterogeneous initial stress from an external file

Usability
4

Just for repetition, the initial state will be included in the prj-file by lines like

251         <parameter>
252             <name>initial_p</name>
253             <type>MeshNode</type>
254             <field_name>pressure</field_name>
255             <mesh>results_equilibrium_ts_10_t_315360000000.000000</mesh>
256         </parameter>
...
311             <initial_condition>initial_p</initial_condition>

for details about stress I would refer to the code documentation (search in the examples of MeshNode for sigma or stress).

Anyway your questions seems to be how to set up a vtu-file (domain mesh) with the initial values, which obviously are not results from a previous OGS run.
So my recommendation is to add point- or cell-data to this mesh by python code like that

import numpy as np      # for operations and data types
from scipy.interpolate import RegularGridInterpolator   # for 3D interpolation
import meshio   # read/write meshes

inputfile="mwe_input.vtu"   # generateStructuredMesh -e hex --lx 10 --ly 10 --lz 10 --nx 20 --ny 20 --nz 20 -o mwe_input.vtu
outputfile="mwe_output.vtu"
    
nnn = np.array([25, 25, 25])  # divisions of data mesh (to be appended to input mesh) in x-, y-, z-direction
xyz = np.array([0,10, 0,10, 0,10])    # start end (x0,x1, y0,y1, z0,z1)
my_data = np.random.standard_normal(size=(nnn[0], nnn[1], nnn[2]))   # generate some data

# data in domain: [X1R,X2R] [Y1R,Y2R] [ZBAS,Z2R] 
NXW = nnn[0]
NYW = nnn[1]
NZW = nnn[2]
x = np.linspace(0.0, xyz[1]-xyz[0], NXW) # shifted to x=0 as mesh does
y = np.linspace(0.0, xyz[3]-xyz[2], NYW) # shifted to y=0 as mesh does
z = np.linspace(0.0, xyz[5]-xyz[4], NZW) # shifted to z=0 as mesh does
print("Data grid nxw, nyw, nzw: ", NXW, NYW, NZW)
print("Data grid x1,x2, y1,y2, zb,zr: ", xyz)

# read input mesh
mesh = meshio.read(inputfile) 
print(f'Input mesh: {len(mesh.points)} points')	# array
print(f'Input mesh: {len(mesh.cells[0].data)} cells')   # dictionary
#print(mesh)

# interpolated data to input mesh ([0,LX]  [0,LY] [0,LZ])
interpolating_function = RegularGridInterpolator((x, y, z), 
                                                 my_data[0:NXW, 0:NYW, 0:NZW], 
                                                 method='nearest',  
                                                 bounds_error=False, 
                                                 fill_value=None)   # method='nearest', 'linear' 

# TODO there may be multiple data on multiple cell blocks
cells_count = len(mesh.cells[0].data)   # MeshIO < 5.0   cells_count = len(mesh.cells[0][1])
centerpoints=np.zeros((cells_count, 3))   # evaluate interpolation at center points
for cell_index, cellblock_cell in enumerate(mesh.cells[0].data):    # MeshIO < 5.0   for cell_index, cellblock_cell in enumerate(mesh.cells[0][1]):
    centerpoints[cell_index] = np.sum(mesh.points[cellblock_cell], axis=0) / len(cellblock_cell)
    
my_data_interpolation = interpolating_function(centerpoints)  

my_data_min = np.min(my_data_interpolation) 
my_data_max = np.max(my_data_interpolation) 
print("Data appended to input mesh are within range {}...{}".format(my_data_min, my_data_max))

mesh.cell_data['my_data'] = [np.array(my_data_interpolation)]   # append data

# write to file
meshio.write(outputfile, mesh)

Without the interpolation stuff the script may be easier, but I do not know what kind of data you have and where they come from.

1 Like