# -*- 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.
"""
def __init__(self, mesh_name, results_dir, pf_type, file_format):
"""
This file is deprecated and would need an overhaul if reactivated in
future versions. For most information, we refer to 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 !!!!!!!!!!!!!!!!!!!!!!! # # #
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)
A_0_r_cg = df.Function(self.V_cg)
A_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)
try:
s1, s2 = self.convertCoordinates(self.mesh.coordinates(), xyz.T)
except IndexError:
print('Primary fields for p2 polynomials are not '
'correct at the moment, aborting ... ')
raise SystemExit
source.config.rho = 1. / np.array(self.sigma)
source.config.d = self.thick
source.config.f = self.f
source.config.current = self.current
source.setLoopMesh(xyz[:, s1])
# lines commented, used for debugging
# np.save('blubb_pos.npy', source.pos)
# np.save('blubb_phi.npy', source.phi)
# np.save('blubb_d.npy', source.ds)
# source.save('blubb')
# raise SystemExit
# lines commented, used for debugging
# ################### calc E-Field ##################################
source.config.ftype = 'E'
E_field = source.calculate(interpolate=False, verbose=False,
maxCPUCount=self.n_procs)
E_0 = E_field[:, s2].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())
A_0_r_cg.vector()[:] = (1.0 / self.omega) * E_0_i_cg.vector()
A_0_i_cg.vector()[:] = -(1.0 / self.omega) * E_0_r_cg.vector()
# ################### calc H-Field ##################################
source.config.ftype = 'H'
H_field = source.calculate(interpolate=False, verbose=False,
maxCPUCount=self.n_procs)
H_0 = H_field[:, s2].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(A_0_r_cg, self.pf_name + '_A_0_real_cg.h5')
self.write_h5(A_0_i_cg, self.pf_name + '_A_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 + '_A_0_real_cg.xml') << A_0_r_cg
df.File(self.pf_name + '_A_0_imag_cg.xml') << A_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 + '_A_0_real_cg.pvd') << A_0_r_cg
df.File(self.pf_name + '_A_0_imag_cg.pvd') << A_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()