Source code for petgem.postprocessing

#!/usr/bin/env python3
# Author:  Octavio Castillo Reyes
# Contact: octavio.castillo@bsc.es
"""Define data postprocessing operations for **PETGEM**."""


# ---------------------------------------------------------------
# Load python modules
# ---------------------------------------------------------------
import numpy as np
import h5py
import meshio
import platform
import shutil
import os
from datetime import datetime
from scipy.spatial import Delaunay
from petsc4py import PETSc

# ---------------------------------------------------------------
# Load petgem modules (BSC)
# ---------------------------------------------------------------
from .common import Print, Timers, measure_all_class_methods
from .parallel import MPIEnvironment, readPetscVector, readPetscMatrix
from .mesh import computeEdges, computeFacesEdges
from .mesh import computeFaces
from .hvfem import computeConnectivityDOFS
from .hvfem import computeJacobian, computeElementOrientation, tetrahedronXYZToXiEtaZeta
from .hvfem import computeBasisFunctions

# ###############################################################
# ################     CLASSES DEFINITION      ##################
# ###############################################################
@measure_all_class_methods
class Postprocessing():
    """Class for postprocessing."""

    def __init__(self):
        """Initialization of a postprocessing class."""
        return

    def run(self, inputSetup):
        """Run a postprocessing task.

        :param obj inputSetup: inputSetup object.
        :return: None
        """
        # ---------------------------------------------------------------
        # Obtain the MPI environment
        # ---------------------------------------------------------------
        parEnv = MPIEnvironment()

        # ---------------------------------------------------------------
        # Initialization
        # ---------------------------------------------------------------
        # Start timer
        Timers()["Postprocessing"].start()

        # Parameters shortcut (for code legibility)
        model = inputSetup.model
        run   = inputSetup.run
        output = inputSetup.output
        out_dir_scratch = output.get('directory_scratch')
        out_dir = output.get('directory')

        # ---------------------------------------------------------------
        # Postprocessing (sequential task)
        # ---------------------------------------------------------------
        if( parEnv.rank == 0 ):
            # ---------------------------------------------------------------
            # Define constants
            # ---------------------------------------------------------------
            num_nodes_per_element = 4
            #num_edges_per_element = 6
            #num_faces_per_element = 4
            #num_nodes_per_face    = 3
            num_edges_per_face    = 3
            #num_nodes_per_edge    = 2
            #num_dimensions        = 3
            basis_order           = run.get('nord')
            num_polarizations     = run.get('num_polarizations')
            num_dof_in_element    = np.int(basis_order*(basis_order+2)*(basis_order+3)/2)
            if (model.get('mode') == 'csem'):
                mode = 'csem'
                data_model = model.get(mode)        # Get data model
                frequency         = data_model.get('source').get('frequency')
            elif (model.get('mode') == 'mt'):
                mode = 'mt'
                data_model = model.get(mode)        # Get data model
                frequency         = data_model.get('frequency')
            omega                 = frequency*2.*np.pi
            mu                    = 4.*np.pi*1e-7
            Const                 = np.sqrt(-1. + 0.j)*omega*mu
            # Get version code
            version_file = open('petgem/VERSION')
            code_version = version_file.read().strip()
            version_file.close()

            # ---------------------------------------------------------------
            # Import mesh file
            # ---------------------------------------------------------------
            mesh_file = model.get('mesh')
            # Import mesh
            mesh = meshio.read(mesh_file)
            # Number of elements
            size = mesh.cells[0][1][:].shape
            nElems = size[0]
            # Number of nodes
            size = mesh.points.shape
            nNodes = size[0]

            # ---------------------------------------------------------------
            # Build data structures for edge elements
            # ---------------------------------------------------------------
            # Compute edges
            elemsE, edgesN = computeEdges(mesh.cells[0][1][:], nElems)
            #nEdges = edgesN.shape[0]

            # Compute faces
            elemsF, facesN = computeFaces(mesh.cells[0][1][:], nElems)
            nFaces = facesN.shape[0]

            # Compute faces-edges
            facesE = computeFacesEdges(elemsF, elemsE, nFaces, nElems)

            # Compute degrees of freedom connectivity
            dofs, dof_edges, dof_faces, _, total_num_dofs = computeConnectivityDOFS(elemsE,elemsF,basis_order)

            # ---------------------------------------------------------------
            # Read receivers file
            # ---------------------------------------------------------------
            # Open receivers_file
            receivers_file = model.get('receivers')
            fileID = h5py.File(receivers_file, 'r')

            # Read receivers
            receivers = fileID.get('data')[()]

            # Number of receivers
            if receivers.ndim == 1:
                nReceivers = 1
            else:
                dim = receivers.shape
                nReceivers = dim[0]

            # ---------------------------------------------------------------
            # Read vector solution
            # ---------------------------------------------------------------
            x = []

            for i in np.arange(num_polarizations):
               input_file = out_dir_scratch + '/x' + str(i) + '.dat'
               x.append(readPetscVector(input_file, communicator=PETSc.COMM_SELF))

            # ---------------------------------------------------------------
            # Compute solution for each point and each polarization mode
            # ---------------------------------------------------------------
            # Allocate array to save field components (Ex, Ey, Ez, Hx, Hy, Hz)
            fields_receivers = []

            for i in np.arange(num_polarizations):
               fields_receivers.append(fieldInterpolator(x[i], mesh.points, mesh.cells[0][1][:], elemsE, edgesN, elemsF, facesE, dofs, receivers, inputSetup))

            # ---------------------------------------------------------------
            # Compute apparent resistivity
            # ---------------------------------------------------------------
            if (mode == 'mt'):
                apparent_res, phase, tipper, impedance = computeImpedance(fields_receivers, omega, mu)

            # ---------------------------------------------------------------
            # Compute solution for entire mesh
            # ---------------------------------------------------------------
            if (output.get('vtk')):
                # Allocate array for points
                scaledCoord = np.zeros((4,3), dtype=np.float)
                local_contribution = np.zeros((num_nodes_per_element,6), dtype=np.complex)
                # Scale factor for tetrahedral element
                scaleFactor = .9999
                # Allocate array to save field components (Ex, Ey, Ez, Hx, Hy, Hz)
                fields_nodes = []

                for i in np.arange(num_polarizations):
                    # Allocate arrays for electromagnetic responses computation
                    fields_tmp = np.zeros((nNodes,6), dtype=np.complex)

                    for j in np.arange(nElems):
                        # Get dofs of element
                        dofsEle = dofs[j, :]

                        # Get indexes of nodes for iEle
                        nodesEle = mesh.cells[0][1][j,:]
                        # Get nodes coordinates for iEle
                        coordEle = mesh.points[nodesEle, :]
                        # Get indexes of faces for iEle
                        facesEle = elemsF[j, :]
                        # Get edges indexes for faces
                        edgesFace = facesE[facesEle, :]
                        # Get indexes of edges for iEle
                        edgesEle = elemsE[j, :]
                        # Get node indexes for edges in i and insert
                        edgesNodesEle = edgesN[edgesEle, :]

                        # Compute jacobian for i
                        jacobian, invjacobian = computeJacobian(coordEle)

                        # Compute global orientation for i
                        edge_orientation, face_orientation = computeElementOrientation(edgesEle,nodesEle,edgesNodesEle,edgesFace)

                        # Compute element centroid
                        centroid = np.sum(coordEle, axis=0)/4.

                        # Apply small scaling to element coordinates
                        scaledCoord[:,0] = centroid[0] + (coordEle[:,0]-centroid[0])*scaleFactor
                        scaledCoord[:,1] = centroid[1] + (coordEle[:,1]-centroid[1])*scaleFactor
                        scaledCoord[:,2] = centroid[2] + (coordEle[:,2]-centroid[2])*scaleFactor

                        # Transform xyz source position to XiEtaZeta coordinates (reference tetrahedral element)
                        XiEtaZeta = tetrahedronXYZToXiEtaZeta(coordEle, scaledCoord)

                        # Compute basis for i
                        basis, curl_basis = computeBasisFunctions(edge_orientation, face_orientation, jacobian, invjacobian, basis_order, XiEtaZeta)

                        # Get global dofs from x vector
                        realField = np.real(x[i].getValues(dofsEle.astype(PETSc.IntType)))
                        imagField = np.imag(x[i].getValues(dofsEle.astype(PETSc.IntType)))

                        # Interpolate field for nodes in i
                        for k in np.arange(num_nodes_per_element):
                            for l in np.arange(num_dof_in_element):
                                # Exyz[i] = Exyz[i] + real_part*basis + imag_part*basis
                                local_contribution[k, 0] += realField[l]*basis[0,l,k] + np.sqrt(-1. + 0.j)*imagField[l]*basis[0,l,k]
                                local_contribution[k, 1] += realField[l]*basis[1,l,k] + np.sqrt(-1. + 0.j)*imagField[l]*basis[1,l,k]
                                local_contribution[k, 2] += realField[l]*basis[2,l,k] + np.sqrt(-1. + 0.j)*imagField[l]*basis[2,l,k]
                                # Hxyz[i] = Hxyz[i] + real_part*curl_basis + imag_part*curl_basis
                                local_contribution[k, 3] += realField[l]*curl_basis[0,l,k] + np.sqrt(-1. + 0.j)*imagField[l]*curl_basis[0,l,k]
                                local_contribution[k, 4] += realField[l]*curl_basis[1,l,k] + np.sqrt(-1. + 0.j)*imagField[l]*curl_basis[1,l,k]
                                local_contribution[k, 5] += realField[l]*curl_basis[2,l,k] + np.sqrt(-1. + 0.j)*imagField[l]*curl_basis[2,l,k]

                        # Add local contribution to global array
                        fields_tmp[nodesEle,:] += local_contribution

                        # Clean variable for next iteration
                        local_contribution[:] = 0.

                    fields_nodes.append(fields_tmp)

                # Compute number of elements that share each node
                num_elems_shared_node = np.zeros(nNodes, dtype=np.float)

                for i in np.arange(nNodes):
                    num_elems_shared_node[i] = np.count_nonzero(mesh.cells[0][1][:] == i)

                # Divide the field of each node by the number of elements that share it
                num_components = 6      # Six electromagnetic field components (Ex, Ey, Ez, Hx, Hy, Hz)
                for i in np.arange(num_polarizations):
                    for j in np.arange(num_components):
                        fields_nodes[i][:,j] /= num_elems_shared_node
                    # Following Maxwell equations, apply constant factor to magnetic field
                    for j in np.arange(3, 6):
                        fields_nodes[i][:,j] /= Const

                # Conductivity model to vtk
                input_file = out_dir_scratch + '/conductivityModel.dat'
                tmpSigmaModel = readPetscMatrix(input_file, communicator=PETSc.COMM_SELF)

                sigmaModel = np.zeros((nElems, 2), dtype=np.float)

                for i in np.arange(nElems):
                    sigmaEle = tmpSigmaModel.getRow(i)[1].real
                    sigmaModel[i:] = sigmaEle

                # Set conductivity model to elements
                mesh.cell_data = {"sigma_horizontal": sigmaModel[:,0], "sigma_vertical": sigmaModel[:,0]}

                # For each polarization and each component, write data to vtk file
                # meshio doesn't support complex number, then fields are decoupled
                # in real and imaginary part
                if (mode == 'csem'):
                    mesh.point_data = {'csem_Ex_real': np.real(fields_nodes[0][:,0]),
                                       'csem_Ey_real': np.real(fields_nodes[0][:,1]),
                                       'csem_Ez_real': np.real(fields_nodes[0][:,2]),
                                       'csem_Hx_real': np.real(fields_nodes[0][:,3]),
                                       'csem_Hy_real': np.real(fields_nodes[0][:,4]),
                                       'csem_Hz_real': np.real(fields_nodes[0][:,5]),
                                       'csem_Ex_imag': np.imag(fields_nodes[0][:,0]),
                                       'csem_Ey_imag': np.imag(fields_nodes[0][:,1]),
                                       'csem_Ez_imag': np.imag(fields_nodes[0][:,2]),
                                       'csem_Hx_imag': np.imag(fields_nodes[0][:,3]),
                                       'csem_Hy_imag': np.imag(fields_nodes[0][:,4]),
                                       'csem_Hz_imag': np.imag(fields_nodes[0][:,5])}
                elif (mode == 'mt'):
                    # Get polarization mode
                    polatization_mode = data_model.get('polarization')
                    if (num_polarizations == 1):
                        mesh.point_data = {'mt_Ex_real_mode'+ polatization_mode: np.real(fields_nodes[0][:,0]),
                                           'mt_Ey_real_mode'+ polatization_mode: np.real(fields_nodes[0][:,1]),
                                           'mt_Ez_real_mode'+ polatization_mode: np.real(fields_nodes[0][:,2]),
                                           'mt_Hx_real_mode'+ polatization_mode: np.real(fields_nodes[0][:,3]),
                                           'mt_Hy_real_mode'+ polatization_mode: np.real(fields_nodes[0][:,4]),
                                           'mt_Hz_real_mode'+ polatization_mode: np.real(fields_nodes[0][:,5]),
                                           'mt_Ex_imag_mode'+ polatization_mode: np.imag(fields_nodes[0][:,0]),
                                           'mt_Ey_imag_mode'+ polatization_mode: np.imag(fields_nodes[0][:,1]),
                                           'mt_Ez_imag_mode'+ polatization_mode: np.imag(fields_nodes[0][:,2]),
                                           'mt_Hx_imag_mode'+ polatization_mode: np.imag(fields_nodes[0][:,3]),
                                           'mt_Hy_imag_mode'+ polatization_mode: np.imag(fields_nodes[0][:,4]),
                                           'mt_Hz_imag_mode'+ polatization_mode: np.imag(fields_nodes[0][:,5])}

                    elif (num_polarizations == 2):
                        mesh.point_data = {'mt_Ex_real_mode'+ polatization_mode[0]: np.real(fields_nodes[0][:,0]),
                                           'mt_Ey_real_mode'+ polatization_mode[0]: np.real(fields_nodes[0][:,1]),
                                           'mt_Ez_real_mode'+ polatization_mode[0]: np.real(fields_nodes[0][:,2]),
                                           'mt_Hx_real_mode'+ polatization_mode[0]: np.real(fields_nodes[0][:,3]),
                                           'mt_Hy_real_mode'+ polatization_mode[0]: np.real(fields_nodes[0][:,4]),
                                           'mt_Hz_real_mode'+ polatization_mode[0]: np.real(fields_nodes[0][:,5]),
                                           'mt_Ex_real_mode'+ polatization_mode[1]: np.real(fields_nodes[1][:,0]),
                                           'mt_Ey_real_mode'+ polatization_mode[1]: np.real(fields_nodes[1][:,1]),
                                           'mt_Ez_real_mode'+ polatization_mode[1]: np.real(fields_nodes[1][:,2]),
                                           'mt_Hx_real_mode'+ polatization_mode[1]: np.real(fields_nodes[1][:,3]),
                                           'mt_Hy_real_mode'+ polatization_mode[1]: np.real(fields_nodes[1][:,4]),
                                           'mt_Hz_real_mode'+ polatization_mode[1]: np.real(fields_nodes[1][:,5]),
                                           'mt_Ex_imag_mode'+ polatization_mode[0]: np.imag(fields_nodes[0][:,0]),
                                           'mt_Ey_imag_mode'+ polatization_mode[0]: np.imag(fields_nodes[0][:,1]),
                                           'mt_Ez_imag_mode'+ polatization_mode[0]: np.imag(fields_nodes[0][:,2]),
                                           'mt_Hx_imag_mode'+ polatization_mode[0]: np.imag(fields_nodes[0][:,3]),
                                           'mt_Hy_imag_mode'+ polatization_mode[0]: np.imag(fields_nodes[0][:,4]),
                                           'mt_Hz_imag_mode'+ polatization_mode[0]: np.imag(fields_nodes[0][:,5]),
                                           'mt_Ex_imag_mode'+ polatization_mode[1]: np.imag(fields_nodes[1][:,0]),
                                           'mt_Ey_imag_mode'+ polatization_mode[1]: np.imag(fields_nodes[1][:,1]),
                                           'mt_Ez_imag_mode'+ polatization_mode[1]: np.imag(fields_nodes[1][:,2]),
                                           'mt_Hx_imag_mode'+ polatization_mode[1]: np.imag(fields_nodes[1][:,3]),
                                           'mt_Hy_imag_mode'+ polatization_mode[1]: np.imag(fields_nodes[1][:,4]),
                                           'mt_Hz_imag_mode'+ polatization_mode[1]: np.imag(fields_nodes[1][:,5])}

                # Write vtk file
                vtk_filename = out_dir + '/' + mode + '_petgemV' + code_version + '_' + str(datetime.today().strftime('%Y-%m-%d')) + '.vtk'
                meshio.write_points_cells(vtk_filename, mesh.points, mesh.cells, mesh.point_data, mesh.cell_data)

            # ---------------------------------------------------------------
            # Save results and data provedance
            # ---------------------------------------------------------------
            # Create output file
            output_file_name = '/' + mode + '_' + 'petgemV' + code_version + '_' + str(datetime.today().strftime('%Y-%m-%d')) + '.h5'
            output_file = out_dir + output_file_name
            fileID = h5py.File(output_file, 'w')

            # Create group for data machine
            machine = fileID.create_group('machine')
            # Name
            uname = platform.uname()
            uname = uname.node + '. ' + uname.release + '. ' + uname.processor
            machine.create_dataset('machine', data=uname)
            # Number of cores
            machine.create_dataset('num_processors', data=parEnv.num_proc)
            # PETGEM version
            machine.create_dataset('petgem_version', data=code_version)

            # Create group for data model
            model_dataprovedance = fileID.create_group('model')
            # Modeling date
            model_dataprovedance.create_dataset('date', data=datetime.today().isoformat())
            # Mesh file
            model_dataprovedance.create_dataset('mesh_file', data=model.get('mesh'))
            # Receivers file
            model_dataprovedance.create_dataset('receivers_file', data=model.get('receivers'))
            # Basis order
            model_dataprovedance.create_dataset('nord', data=run.get('nord'))
            # Cuda support
            model_dataprovedance.create_dataset('cuda', data=run.get('cuda'))
            # vtk output
            model_dataprovedance.create_dataset('vtk', data=output.get('vtk'))
            # Modeling mode
            model_dataprovedance.create_dataset('mode', data=mode)
            # Number of polarizations
            model_dataprovedance.create_dataset('num_polarizations', data=num_polarizations)
            # Create data model
            if (mode == 'csem'):
                if (run.get('conductivity_from_file')):
                    model_dataprovedance.create_dataset('sigma', data=data_model.get('sigma').get('file'))
                else:
                    model_dataprovedance.create_dataset('sigma_horizontal', data=data_model.get('sigma').get('horizontal'))
                    model_dataprovedance.create_dataset('sigma_vertical', data=data_model.get('sigma').get('vertical'))
                model_dataprovedance.create_dataset('frequency', data=data_model.get('source').get('frequency'))
                model_dataprovedance.create_dataset('source_position', data=np.asarray(data_model.get('source').get('position'), dtype=np.float))
                model_dataprovedance.create_dataset('source_azimuth', data=data_model.get('source').get('azimuth'))
                model_dataprovedance.create_dataset('source_dip', data=data_model.get('source').get('dip'))
                model_dataprovedance.create_dataset('source_current', data=data_model.get('source').get('current'))
                model_dataprovedance.create_dataset('source_length', data=data_model.get('source').get('length'))
            elif (mode == 'mt'):
                if (run.get('conductivity_from_file')):
                    model_dataprovedance.create_dataset('sigma', data=data_model.get('sigma').get('file'))
                else:
                    model_dataprovedance.create_dataset('sigma_horizontal', data=data_model.get('sigma').get('horizontal'))
                    model_dataprovedance.create_dataset('sigma_vertical', data=data_model.get('sigma').get('vertical'))
                model_dataprovedance.create_dataset('frequency', data=data_model.get('frequency'))
                model_dataprovedance.create_dataset('polarization', data=data_model.get('polarization'))

            # Write electromagnetic responses
            if (mode == 'csem'):
                # Electric fields
                E_fields = model_dataprovedance.create_group('E-fields')
                E_fields.create_dataset('x', data=fields_receivers[0][:,0])
                E_fields.create_dataset('y', data=fields_receivers[0][:,1])
                E_fields.create_dataset('z', data=fields_receivers[0][:,2])

                # Magnetic fields
                H_fields = model_dataprovedance.create_group('H-fields')
                H_fields.create_dataset('x', data=fields_receivers[0][:,3])
                H_fields.create_dataset('y', data=fields_receivers[0][:,4])
                H_fields.create_dataset('z', data=fields_receivers[0][:,5])

            elif (mode == 'mt'):
                list_modes = data_model.get('polarization')
                for i in np.arange(num_polarizations):
                    mode_E_i = 'E-fields_mode_' + list_modes[i]
                    mode_H_i = 'H-fields_mode_' + list_modes[i]

                    E_fields = model_dataprovedance.create_group(mode_E_i)
                    E_fields.create_dataset('x', data=fields_receivers[i][:,0])
                    E_fields.create_dataset('y', data=fields_receivers[i][:,1])
                    E_fields.create_dataset('z', data=fields_receivers[i][:,2])

                    # Magnetic fields
                    H_fields = model_dataprovedance.create_group(mode_H_i)
                    H_fields.create_dataset('x', data=fields_receivers[i][:,3])
                    H_fields.create_dataset('y', data=fields_receivers[i][:,4])
                    H_fields.create_dataset('z', data=fields_receivers[i][:,5])

                imp = model_dataprovedance.create_group('impedance')
                imp.create_dataset('xx', data=impedance[0])
                imp.create_dataset('xy', data=impedance[1])
                imp.create_dataset('yx', data=impedance[2])
                imp.create_dataset('yy', data=impedance[3])

                app_res = model_dataprovedance.create_group('apparent_resistivity')
                app_res.create_dataset('xx', data=apparent_res[0])
                app_res.create_dataset('xy', data=apparent_res[1])
                app_res.create_dataset('yx', data=apparent_res[2])
                app_res.create_dataset('yy', data=apparent_res[3])

                pha = model_dataprovedance.create_group('phase')
                pha.create_dataset('xx', data=phase[0])
                pha.create_dataset('xy', data=phase[1])
                pha.create_dataset('yx', data=phase[2])
                pha.create_dataset('yy', data=phase[3])

                tip = model_dataprovedance.create_group('tipper')
                tip.create_dataset('x', data=tipper[0])
                tip.create_dataset('y', data=tipper[1])

            # Close file
            fileID.close()

            # Remove temporal directory
            if (output.get('remove_scratch')):
                shutil.rmtree(out_dir_scratch)
            else:
                files_in_directory = os.listdir(out_dir_scratch)
                filtered_files = [file for file in files_in_directory if (file.endswith(".dat") or file.endswith(".info"))]
                for file in filtered_files:
                    path_to_file = os.path.join(out_dir_scratch, file)
                    os.remove(path_to_file)

        # Stop timer
        Timers()["Postprocessing"].stop()

        return


[docs]def fieldInterpolator(solution_vector, nodes, elemsN, elemsE, edgesN, elemsF, facesE, dof_connectivity, points, inputSetup): """Interpolate electromagnetic field for a set of 3D points. :param ndarray-petsc solution_vector: vector field to be interpolated :param ndarray nodes: nodal coordinates :param ndarray elemsN: elements-node connectivity with dimensions = (number_elements, 4) :param ndarray elemsE: elements-edge connectivity with dimensions = (number_elements, 6) :param ndarray edgesN: edge-node connectivity with dimensions = (number_edges, 2) :param ndarray elemsF: element-faces connectivity with dimensions = (number_elements, 4) :param ndarray facesE: face-edges connectivity with dimensions = (number_faces, 3) :param ndarray dof_connectivity: local/global dofs list for elements, dofs index on edges, dofs index on faces, dofs index on volumes, total number of dofs :param ndarray points: point coordinates :param obj inputSetup: inputSetup object. :return: electromagnetic fields for a set of 3D points :rtype: ndarray and int """ # --------------------------------------------------------------- # Initialization # --------------------------------------------------------------- # Parameters shortcut (for code legibility) model = inputSetup.model run = inputSetup.run basis_order = run.get('nord') if (model.get('mode') == 'csem'): mode = 'csem' data_model = model.get(mode) # Get data model frequency = data_model.get('source').get('frequency') elif (model.get('mode') == 'mt'): mode = 'mt' data_model = model.get(mode) # Get data model frequency = data_model.get('frequency') omega = frequency*2.*np.pi mu = 4.*np.pi*1e-7 Const = np.sqrt(-1. + 0.j)*omega*mu # Number of elements size = elemsN.shape nElems = size[0] # Number of nodes #size = nodes.shape #nNodes = size[0] # Num dof per element num_dof_in_element = np.int(basis_order*(basis_order+2)*(basis_order+3)/2) # Number of points if points.ndim == 1: nPoints = 1 else: dim = points.shape nPoints = dim[0] # Find where receivers are located tri = Delaunay(nodes) # Overwrite Delaunay structure with mesh_file connectivity and points tri.simplices = elemsN tri.vertices = elemsN # Find out which tetrahedral element points are in points_in_elements = tri.find_simplex(points, bruteforce=True, tol=1.e-12) # Determine if all points were found idx = np.where(np.logical_or(points_in_elements>nElems, points_in_elements<0))[0] # If idx is not empty, there are points outside the domain if idx.size != 0: Print.master(' The following receivers were not located and will not be taken into account ' + str(idx)) # Update number of receivers nPoints = nPoints - len(idx) if nPoints == 0: Print.master(' No point has been found. Nothing to do. Aborting') exit(-1) # Remove idx from points matrix points = np.delete(points, idx, axis=0) # Remove idx from points_in_elements points_in_elements = np.delete(points_in_elements, idx, axis=0) indx_ele = points_in_elements # Allocate array fields = np.zeros((nPoints, 6), dtype=np.complex) # Interpolate electromagnetic field for all points for i in np.arange(nPoints): # Get element index iEle = indx_ele[i] # Get dofs of element container dofsEle = dof_connectivity[iEle, :] # Get receiver coordinates coordPoints = points[i,:] # Get indexes of nodes for iEle nodesEle = elemsN[iEle,:] # Get nodes coordinates for iEle coordEle = nodes[nodesEle, :] # Get indexes of faces for iEle facesEle = elemsF[iEle, :] # Get edges indexes for faces edgesFace = facesE[facesEle, :] # Get indexes of edges for iEle edgesEle = elemsE[iEle, :] # Get node indexes for edges in i and insert edgesNodesEle = edgesN[edgesEle, :] # Compute jacobian for iEle jacobian, invjacobian = computeJacobian(coordEle) # Compute global orientation for iEle edge_orientation, face_orientation = computeElementOrientation(edgesEle,nodesEle,edgesNodesEle,edgesFace) # Transform xyz source position to XiEtaZeta coordinates (reference tetrahedral element) XiEtaZeta = tetrahedronXYZToXiEtaZeta(coordEle, coordPoints) # Compute basis for i basis, curl_basis = computeBasisFunctions(edge_orientation, face_orientation, jacobian, invjacobian, basis_order, XiEtaZeta) # Get global dofs from x vector realField = np.real(solution_vector.getValues(dofsEle.astype(PETSc.IntType))) imagField = np.imag(solution_vector.getValues(dofsEle.astype(PETSc.IntType))) for j in np.arange(num_dof_in_element): # Exyz[k] = Exyz[k] + real_part*basis + imag_part*basis fields[i, 0] += realField[j]*basis[0,j] + np.sqrt(-1. + 0.j)*imagField[j]*basis[0,j] fields[i, 1] += realField[j]*basis[1,j] + np.sqrt(-1. + 0.j)*imagField[j]*basis[1,j] fields[i, 2] += realField[j]*basis[2,j] + np.sqrt(-1. + 0.j)*imagField[j]*basis[2,j] # Hxyz[k] = Hxyz[k] + real_part*curl_basis + imag_part*curl_basis fields[i, 3] += realField[j]*curl_basis[0,j] + np.sqrt(-1. + 0.j)*imagField[j]*curl_basis[0,j] fields[i, 4] += realField[j]*curl_basis[1,j] + np.sqrt(-1. + 0.j)*imagField[j]*curl_basis[1,j] fields[i, 5] += realField[j]*curl_basis[2,j] + np.sqrt(-1. + 0.j)*imagField[j]*curl_basis[2,j] # Following Maxwell equations, apply constant factor to compute magnetic field fields[:,3::] /= Const return fields
[docs]def computeImpedance(fields, omega, mu): """Compute apparent resistiviy, phase, tipper and impedance for MT mode. :param list fields: list of numpy arrays with electromagnetic fields with dimensions = (number_polarizations, number_receivers, number_EM_components) :param float omega: angular frequency :param float mu: medium permeability. :return: apparent resistivity, phase, tipper and impedance. :rtype: ndarray """ # Field components (Ex, Ey, Hx, Hy and Hz components por x and y-directions in polarization) E1x = fields[0][:,0] E1y = fields[0][:,1] E2x = fields[1][:,0] E2y = fields[1][:,1] H1x = fields[0][:,3] H1y = fields[0][:,4] H1z = fields[0][:,5] H2x = fields[1][:,3] H2y = fields[1][:,4] H2z = fields[1][:,5] # Minor determinants m1 = H2y * H2x * H1y - H2y * H2y * H1x m2 = H2x * H2x * H1y - H2y * H2x * H1x m3 = E2x * H2x * H1y - H2y * E2x * H1x m4 = H2y * E2y * H1y - H2y * H2y * E1y m5 = H2x * E2y * H1y - H2y * H2x * E1y m6 = H2y * H2x * E1y - E2y * H2y * H1x m7 = H2x * H2x * E1y - E2y * H2x * H1x # Major determinants d = H1x * m1 - H1y * m2 d1 = E1x * m1 - H1y * m3 d2 = H1x * m3 - E1x * m2 d3 = H1x * m4 - H1y * m5 d4 = H1x * m6 - H1y * m7 # Impedance impedance = [] impedance.append(d1/d) # xx-impedance impedance.append(d2/d) # xy-impedance impedance.append(d3/d) # yx-impedance impedance.append(d4/d) # yy-impedance # Aparent resistivity coef = np.float(1.) / (mu * omega) apparent_resistivity = [] apparent_resistivity.append(coef * np.abs(impedance[0])**2) # xx-apparent resistivity apparent_resistivity.append(coef * np.abs(impedance[1])**2) # xy-apparent resistivity apparent_resistivity.append(coef * np.abs(impedance[2])**2) # yx-apparent resistivity apparent_resistivity.append(coef * np.abs(impedance[3])**2) # yy-apparent resistivity # Phase coef = np.float(180.) / np.pi phase = [] phase.append((coef * np.arctan(-np.imag(impedance[0]) / np.real(impedance[0])))%360) # xx-phase phase.append((coef * np.arctan(-np.imag(impedance[1]) / np.real(impedance[1])))%360) # xy-phase phase.append((coef * np.arctan(-np.imag(impedance[2]) / np.real(impedance[2])))%360) # yx-phase phase.append((coef * np.arctan(-np.imag(impedance[3]) / np.real(impedance[3])))%360) # yy-phase # Tipper --> solve for T = (Tx, Ty) such that Hz = <T, H> with H = (Hx, Hy) coef = H1x * H2y tmp1 = (H2z/H2y - (H1z*H2x)/coef) / (np.float(1.) - (H1y*H2x)/coef) tmp2 = (H1z - tmp1 * H1y) / H1x tipper = [] tipper.append(tmp2) # x-tipper tipper.append(tmp1) # y-tipper return apparent_resistivity, phase, tipper, impedance
# ############################################################### # ################ FUNCTIONS DEFINITION ################# # ###############################################################
[docs]def unitary_test(): """Unitary test for postprocessing.py script."""
# ############################################################### # ################ MAIN ################# # ############################################################### if __name__ == '__main__': # --------------------------------------------------------------- # Run unitary test # --------------------------------------------------------------- unitary_test()