diff --git a/gprMax/geometry_outputs.py b/gprMax/geometry_outputs.py index f36b2d44..a4dbcb5e 100644 --- a/gprMax/geometry_outputs.py +++ b/gprMax/geometry_outputs.py @@ -23,10 +23,13 @@ except ImportError: import logging from pathlib import Path from struct import pack +from xml.etree.ElementTree import indent import gprMax.config as config import h5py import numpy as np +from evtk.hl import rectilinearToVTK +from evtk.hl import linesToVTK from ._version import __version__ from .cython.geometry_outputs import (define_fine_geometry, @@ -124,6 +127,83 @@ class GeometryView: self.filename = Path(*parts[:-1], self.filename + config.get_model_config().appendmodelnumber) self.filename = self.filename.with_suffix(self.fileext) + + def evtk_line(self): + + G = self.G + # line counter + lc = 0 + # point counter + pc = 0 + + n_x_lines = self.nx * (self.ny + 1) * (self.nz + 1) + n_y_lines = self.ny * (self.nx + 1) * (self.nz + 1) + n_z_lines = self.nz * (self.nx + 1) * (self.ny + 1) + + n_lines = n_x_lines + n_y_lines + n_z_lines + + # a line is defined by 2 points. Each line must have 2 new points + n_points = 2 * n_lines + + x = np.zeros(n_points) + y = np.zeros(n_points) + z = np.zeros(n_points) + + l = np.zeros(n_lines) + + for i in range(self.nx): + for j in range(self.ny): + for k in range(self.nz): + + # x yee cell edge + # line the line property to the relevent material + + l[lc] = G.ID[0][i, j, k] + # set the starting point position of the edge + x[pc], y[pc], z[pc] = i * self.dx, j * self.dy, k * self.dz + # next point + pc += 1 + # set the end point position of the edge + x[pc], y[pc], z[pc] = (i + 1) * self.dx, j * self.dy, k * self.dz + # next point + pc += 1 + + # next line + lc += 1 + + # y yee cell edge + l[lc] = self.G.ID[1, i, j, k] + x[pc], y[pc], z[pc] = i * self.dx, j * self.dy, k * self.dz + pc += 1 + x[pc], y[pc], z[pc] = i * self.dx, (j + 1) * self.dy, k * self.dz + pc += 1 + lc += 1 + + # z yee cell edge + l[lc] = self.G.ID[2, i, j, k] + x[pc], y[pc], z[pc] = i * self.dx, j * self.dy, k * self.dz + pc += 1 + x[pc], y[pc], z[pc] = i * self.dx, j * self.dy, (k + 1) * self.dz + pc += 1 + lc += 1 + + comments = [ "comment 1", "comment 2" ] + linesToVTK('./linesToVTK', x, y, z, cellData = {"Yee materials" : l}) + + def evtk_voxels(self):# have a go with evtk-hl library + G = self.G + lx = G.solid.shape[0] * G.dx + ly = G.solid.shape[1] * G.dy + lz = G.solid.shape[2] * G.dz + + # Coordinates + x = np.arange(0, lx + 0.1 * G.dx, G.dx, dtype='float64') + y = np.arange(0, ly + 0.1 * G.dx, G.dy, dtype='float64') + z = np.arange(0, lz + 0.1 * G.dx, G.dz, dtype='float64') + + fp = './test-evtk-hl' + rectilinearToVTK(fp, x, y, z, cellData = {"solid" : G.solid}, comments = ['hi', 'world']) + def write_vtk(self, G, pbar): """Writes the geometry information to a VTK file. Either ImageData (.vti) for a per-cell geometry view, or @@ -173,26 +253,12 @@ class GeometryView: f.write('\n\n\n\n'.encode('utf-8')) self.write_gprmax_info(f, G) - # have a go with evtk-hl library - from time import perf_counter - t0 = perf_counter() - import os - from evtk.hl import rectilinearToVTK - - lx = G.solid.shape[0] * G.dx - ly = G.solid.shape[1] * G.dy - lz = G.solid.shape[2] * G.dz - - # Coordinates - x = np.arange(0, lx + 0.1 * G.dx, G.dx, dtype='float64') - y = np.arange(0, ly + 0.1 * G.dx, G.dy, dtype='float64') - z = np.arange(0, lz + 0.1 * G.dx, G.dz, dtype='float64') - - fp = './test-evtk-hl' - rectilinearToVTK(fp, x, y, z, cellData = {"solid" : G.solid}, comments = ['hi', 'world']) - print(f'time to write using evtk-hl {perf_counter() - t0}') + self.evtk_voxels() elif self.fileext == '.vtp': + + self.evtk_line() + with open(self.filename, 'wb') as f: f.write('\n'.encode('utf-8')) f.write(f'\n'.encode('utf-8')) diff --git a/tests/other_codes/evtk_test.py b/tests/other_codes/evtk_test.py new file mode 100644 index 00000000..13c06270 --- /dev/null +++ b/tests/other_codes/evtk_test.py @@ -0,0 +1,83 @@ +# Copyright (C) 2015-2021, John Hartley +# +# This module is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. +# To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/. + +from pathlib import Path +import sys +import gprMax +from user_libs.antennas.GSSI import antenna_like_GSSI_400 + +import numpy as np + +# file path step +fn = Path(__file__) +parts = fn.parts + +ratio = 5 +dl = 5e-3 + +# cells +# default number of pml cells +pml_cells = 10 +# distance between model and pml cells +pml_gap = 15 +# number of cells between the Inner Surface and the Outer Surface of the sub-grid +is_os_gap = 4 +# size of the sub-gridded region +sub_gridded_region = 3 +# domain size +extent = 3 + 2 * (pml_cells + pml_gap + 4) + +# domain extent +x = dl * extent +y = x +z = x + +tw = 1e-9 + +scene = gprMax.Scene() + +title_gpr = gprMax.Title(name=fn.name) +dxdydz = gprMax.Discretisation(p1=(dl, dl, dl)) +domain = gprMax.Domain(p1=(x, y, z)) +time_window = gprMax.TimeWindow(time=tw) + +scene.add(domain) +scene.add(title_gpr) +scene.add(dxdydz) +scene.add(time_window) + + +# plastic box in sub grid +material = gprMax.Material(er=3, mr=1, se=0, sm=0, id='plastic') +scene.add(material) + +plastic_box = gprMax.Box(p1=(20*dl, 20*dl, 20*dl), p2=(40*dl, 40*dl, 30*dl), material_id='plastic') +scene.add(plastic_box) + +pec_box = gprMax.Box(p1=(20*dl, 20*dl, 30*dl), p2=(40*dl, 40*dl, 40*dl), material_id='pec') +scene.add(pec_box) + +pec_plate = gprMax.Plate(p1=(20*dl, 20*dl, 20*dl), p2=(40*dl, 40*dl, 20*dl), material_id='pec') +scene.add(pec_plate) + +# create a geometry view of the main grid and the sub grid stitched together +gv1 = gprMax.GeometryView(p1=(0, 0, 0), + p2=domain.props.p1, + dl=dl, + filename=fn.with_suffix('').parts[-1], + output_type='f') + +# create a geometry view of the main grid and the sub grid stitched together +gv2 = gprMax.GeometryView(p1=(0, 0, 0), + p2=domain.props.p1, + dl=dl, + filename=fn.with_suffix('').parts[-1], + output_type='n') + + +scene.add(gv1) +scene.add(gv2) + +gprMax.run(scenes=[scene], n=1, geometry_only=False, outputfile=fn)