# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
import dolfin as df
import os
import numpy as np
import sys
import time
import getopt
from comet import pyhed as ph
from mpi4py import MPI
"""
This file is always called in serial from the root process of an mpi run,
otherwise a crash would occur due to a missing 'MPI-Parent' for synchronization
"""
[docs]
class Pyhed_calc:
"""
Primary field calculation class which starts calculating primary fields
with *pyhed*-internal multiprocessing from a serial process to ensure
correct dof ordering.
"""
def __init__(self, mesh_name, results_dir, pf_type, file_format):
"""
This module needs an update in future versions. There are no specific
docs for this file, but most content can be understood by reading the
docs of the *PrimaryFields* class.
"""
self.mpi_cw = df.MPI.comm_world
self.mesh_name = mesh_name
self.r_dir = results_dir
self.pf_type = pf_type
self.file_format = file_format
self.read_serial_calculation_parameters()
if sys.version_info[:2] < (3, 0):
self.pf_name = self.pf_name.encode('ascii', 'ignore')
self.pvd_flag = False
self.calc_fields()
[docs]
def read_serial_calculation_parameters(self):
import json
path = os.getcwd()
os.chdir(self.r_dir + '/primary_fields/' + self.pf_type)
with open('pf_dict.json', 'r') as thefile:
self.__dict__.update(json.load(thefile))
if self.tx == 'import':
self.tx = np.load('source_path.npy')
os.chdir(path)
[docs]
def write_h5(self, var, fname):
f = df.HDF5File(self.mpi_cw, fname, "w")
f.write(var, "/data")
f.close()
[docs]
def read_h5(self, var, fname):
f = df.HDF5File(self.mpi_cw, fname, "r")
f.read(var, "/data")
f.close()
[docs]
def convertCoordinates(self, gimli, dolfin):
"""
Input: arr1, arr2: two coordinate lists of same shape (n, 3)
which contains the same coordinates but in a diffrent order.
output: arr1_arr2, arr2_arr1: index arrays which converts coordinates
from input1 to input2 and from input2 to input1.
"""
def indexAndSorting(c_):
posview = c_.ravel().view(
np.dtype((np.void, c_.dtype.itemsize*c_.shape[1])))
_, unique_idx = np.unique(posview, return_index=True)
return unique_idx, unique_idx.argsort()
pg_idx, pg_sort = indexAndSorting(gimli)
df_idx, df_sort = indexAndSorting(dolfin)
return pg_idx[df_sort], df_idx[pg_sort]
[docs]
def calc_fields(self):
path = os.getcwd()
self.mesh = df.Mesh()
# # # !!!!!!!!!!!!!!!!!!!!!!!!!!!! HACKED !!!!!!!!!!!!!!!!!!!!!!! # # #
if os.path.isfile('../meshes/_h5/' + self.mesh_name + '.h5'):
hdf5 = df.HDF5File(self.mpi_cw,
'../meshes/_h5/' + self.mesh_name + '.h5', "r")
elif os.path.isfile('../../meshes/_h5/' + self.mesh_name + '.h5'):
hdf5 = df.HDF5File(self.mpi_cw, '../../meshes/_h5/' +
self.mesh_name + '.h5', "r")
elif os.path.isfile('meshes/_h5/' + self.mesh_name + '.h5'):
hdf5 = df.HDF5File(self.mpi_cw, 'meshes/_h5/' +
self.mesh_name + '.h5', "r")
hdf5.read(self.mesh, '/mesh', False)
# # # !!!!!!!!!!!!!!!!!!!!!!!!!!!! HACKED !!!!!!!!!!!!!!!!!!!!!!! # # #
if self.pf_p == 2:
print(' - using 2nd polynomial order for primary fields - ')
self.V_cg = df.VectorFunctionSpace(self.mesh, 'CG', self.pf_p)
t1 = time.time()
E_0_r_cg = df.Function(self.V_cg)
E_0_i_cg = df.Function(self.V_cg)
H_0_r_cg = df.Function(self.V_cg)
H_0_i_cg = df.Function(self.V_cg)
xyz = (self.V_cg.tabulate_dof_coordinates().reshape(-1, 3)[0::3, :]).T
if self.s_type == 'hed':
source = ph.loop.buildDipole(self.origin, length=self.length,
angle=np.deg2rad(self.azimuth))
print('########## Primary "' + self.pf_type + '_' + self.s_type +
'" field is being calculated ... ###########\nHED origin:')
print(' Origin' + ' --> ', self.origin)
elif self.s_type == 'line':
source = ph.loop.buildLine(self.start, self.stop,
num_segs=self.n_segs)
print('########## Primary "' + self.pf_type + '_' + self.s_type +
'" field is being calculated ... ###########\n'
'Line source parameters:')
print(' Start_point' + ' --> ', self.start)
print(' Stop_point' + ' --> ', self.stop)
print(' No. of segments' + ' --> ', self.n_segs)
elif self.s_type == 'loop_c':
source = ph.loop.buildCircle(self.r, P=self.origin,
num_segs=self.n_segs)
print('########## Primary "' + self.pf_type + '_' + self.s_type +
'" field is being calculated ... ###########\n'
'Circular loop source parameters:')
print(' Origin' + ' --> ', self.origin)
print(' Radius' + ' --> ', self.r)
print(' No. of segments' + ' --> ', self.n_segs)
elif self.s_type == 'loop_r':
if df.near(self.start[0], self.stop[0], eps=1e-4) or \
df.near(self.start[1], self.stop[1], eps=1e-4):
print('Fatal error: Given start and stop coordinates do not '
'match to a rectangular loop source: x or y start or '
'stop coordinates are identical')
raise SystemExit
S_1 = [self.start[0], self.start[1], 0.0]
S_2 = [self.stop[0], self.start[1], 0.0]
S_3 = [self.stop[0], self.stop[1], 0.0]
S_4 = [self.start[0], self.stop[1], 0.0]
source = ph.loop.buildRectangle(np.array((S_1, S_2,
S_3, S_4)),
num_segs=self.n_segs)
print('########## Primary "' + self.pf_type + '_' + self.s_type +
'" field is being calculated ... ###########\n'
'Rectangular loop source parameters:')
print(' Corner 1' + ' --> ', S_1)
print(' Corner 2' + ' --> ', S_2)
print(' Corner 3' + ' --> ', S_3)
print(' Corner 4' + ' --> ', S_4)
print(' No. of segments' + ' --> ', self.n_segs/4)
elif self.s_type == 'path':
source = ph.loop.buildLoop(self.tx[0], num_segs=self.n_segs,
grounded=self.grounding[0])
print('########## Primary "' + self.pf_type + '_' + self.s_type +
'" field is being calculated ... ###########\n'
'Source from path parameters:')
print(' closed loop source' + ' --> ', not self.grounding[0])
else:
print('!!! Fatal Error, supported "s_type" parameters for '
'Halfspace" primary fields are: "HED, "Line", "Loop_C"'
' or "Loop_R"')
raise SystemExit
os.chdir(self.r_dir + '/primary_fields/' + self.pf_type)
selection = np.zeros(xyz.shape[1], dtype=bool)
xlim = [np.min(xyz[0, :]), np.max(xyz[0, :])]
ylim = [np.min(xyz[1, :]), np.max(xyz[1, :])]
zlim = [np.min(xyz[2, :]), np.max(xyz[2, :])]
self.boundary_eps = 1e-2
for j in range(xyz.shape[1]):
if df.near(xyz[0, j], xlim[0], self.boundary_eps) or \
df.near(xyz[0, j], xlim[1], self.boundary_eps):
selection[j] = True
elif (df.near(xyz[1, j], ylim[0], self.boundary_eps) or
df.near(xyz[1, j], ylim[1], self.boundary_eps)):
selection[j] = True
elif (df.near(xyz[2, j], zlim[0], self.boundary_eps) or
df.near(xyz[2, j], zlim[1], self.boundary_eps)):
selection[j] = True
else:
pass
source.setLoopMesh(xyz[:, selection])
source.config.rho = 1. / np.array(self.sigma)
source.config.d = self.thick
source.config.f = self.f
source.config.current = self.current
# ################### calc E-Field ##################################
source.config.ftype = 'E'
E_field = source.calculate(interpolate=False, verbose=False,
maxCPUCount=self.n_procs)
E_0 = np.zeros((xyz.shape[1], 3), dtype=complex)
E_0[selection, :] = E_field[:, :].T
E_0_r_cg.vector()[:] = (np.hstack((E_0[:, 0].real[np.newaxis].T,
np.hstack((E_0[:, 1].real[np.newaxis].T,
E_0[:, 2].real[np.newaxis].T))
)).ravel())
E_0_i_cg.vector()[:] = (np.hstack((E_0[:, 0].imag[np.newaxis].T,
np.hstack((E_0[:, 1].imag[np.newaxis].T,
E_0[:, 2].imag[np.newaxis].T))
)).ravel())
# ################### calc H-Field ##################################
source.config.ftype = 'H'
H_field = source.calculate(interpolate=False, verbose=False,
maxCPUCount=self.n_procs)
H_0 = np.zeros((xyz.shape[1], 3), dtype=complex)
H_0[selection, :] = H_field[:, :].T
H_0_r_cg.vector().set_local((np.hstack((H_0[:, 0].real[
np.newaxis].T,
np.hstack((H_0[:, 1].real[
np.newaxis].T,
H_0[:, 2].real[
np.newaxis].T))
)).ravel()))
H_0_r_cg.vector().apply('insert')
H_0_i_cg.vector().set_local((np.hstack((H_0[:, 0].imag[
np.newaxis].T,
np.hstack((H_0[:, 1].imag[
np.newaxis].T,
H_0[:, 2].imag[
np.newaxis].T))
)).ravel()))
H_0_i_cg.vector().apply('insert')
print('... saving primary field in "primary_fields" directory ...')
if self.file_format == 'h5':
self.write_h5(E_0_r_cg, self.pf_name + '_E_0_real_cg.h5')
self.write_h5(E_0_i_cg, self.pf_name + '_E_0_imag_cg.h5')
self.write_h5(H_0_r_cg, self.pf_name + '_H_0_real_cg.h5')
self.write_h5(H_0_i_cg, self.pf_name + '_H_0_imag_cg.h5')
elif self.file_format == 'xml':
df.File(self.pf_name + '_E_0_real_cg.xml') << E_0_r_cg
df.File(self.pf_name + '_E_0_imag_cg.xml') << E_0_i_cg
df.File(self.pf_name + '_H_0_real_cg.xml') << H_0_r_cg
df.File(self.pf_name + '_H_0_imag_cg.xml') << H_0_i_cg
if self.pvd_flag:
df.File(self.pf_name + '_E_0_real_cg.pvd') << E_0_r_cg
df.File(self.pf_name + '_E_0_imag_cg.pvd') << E_0_i_cg
df.File(self.pf_name + '_H_0_real_cg.pvd') << H_0_r_cg
df.File(self.pf_name + '_H_0_imag_cg.pvd') << H_0_i_cg
# ################### calc H-Field if possible ######################
print('... done.')
t2 = time.time() - t1
print('##############################################################')
print(' Computation time (s) for primary fields : ', t2)
print('##############################################################')
os.chdir(path)
if __name__ == "__main__":
myopts, args = getopt.getopt(sys.argv[1:], "m:d:t:f:")
for oo, aa in myopts:
if oo == '-m':
mesh = str(aa)
elif oo == '-d':
results_dir = str(aa)
elif oo == '-t':
pf_type = str(aa)
elif oo == '-f':
file_format = str(aa)
else:
print("Usage: %s -m mesh_name, -d results_dir"
", -t pf_type -f file_format" % sys.argv[0])
parent = MPI.Comm.Get_parent()
Pyhed_calc(mesh, results_dir, pf_type, file_format)
parent.Barrier()
parent.Disconnect()