Source code for petgem.preprocessing

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


# ---------------------------------------------------------------
# Load python modules
# ---------------------------------------------------------------
import numpy as np
import h5py
import meshio
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, createSequentialDenseMatrixWithArray
from .parallel import writeParallelDenseMatrix, createSequentialVectorWithArray
from .parallel import writePetscVector
from .mesh import computeEdges, computeBoundaryEdges, computeFacesEdges
from .mesh import computeFaces, computeBoundaryFaces
from .mesh import computeBoundaryElements, computeBoundaries, computeFacePlane
from .hvfem import computeConnectivityDOFS


# ###############################################################
# ################     CLASSES DEFINITION      ##################
# ###############################################################
@measure_all_class_methods
class Preprocessing():
    """Class for preprocessing."""

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

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

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

        # Start timer
        Timers()["Preprocessing"].start()

        # ---------------------------------------------------------------
        # Preprocessing (sequential task)
        # ---------------------------------------------------------------
        if( parEnv.rank == 0 ):
            # Parameters shortcut (for code legibility)
            model = inputSetup.model
            run = inputSetup.run
            output = inputSetup.output
            out_dir = output.get('directory_scratch')
            # Compute number of dofs per element
            basis_order = run.get('nord')
            num_dof_in_element = np.int(basis_order*(basis_order+2)*(basis_order+3)/2)
            if (model.get('mode') == 'csem'):
                mode = 'csem'
            elif (model.get('mode') == 'mt'):
                mode = 'mt'
            # Get data model
            data_model = model.get(mode)

            # ---------------------------------------------------------------
            # 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]

            # ---------------------------------------------------------------
            # Preprocessing nodal coordinates
            # ---------------------------------------------------------------
            Print.master('     Nodal coordinates')

            # Build coordinates in PETGEM format where each row
            # represent the xyz coordinates of the 4 tetrahedral element
            num_dimensions = 3
            num_nodes_per_element = 4
            data = mesh.points[mesh.cells[0][1][:], :]
            data = data.reshape(nElems, num_dimensions*num_nodes_per_element)

            # Get matrix dimensions
            size = data.shape

            # Build PETSc structures
            matrix = createSequentialDenseMatrixWithArray(size[0], size[1], data)

            # Build path to save the file
            out_path = out_dir + '/nodes.dat'

            # Write PETGEM nodes in PETSc format
            writeParallelDenseMatrix(out_path, matrix, communicator=PETSc.COMM_SELF)
            # Remove temporal matrix
            del matrix

            # ---------------------------------------------------------------
            # Preprocessing mesh connectivity
            # ---------------------------------------------------------------
            Print.master('     Mesh connectivity')

            # Get matrix dimensions
            size = mesh.cells[0][1][:].shape

            # Build PETSc structures
            matrix = createSequentialDenseMatrixWithArray(size[0], size[1], mesh.cells[0][1][:])

            # Build path to save the file
            out_path = out_dir + '/meshConnectivity.dat'

            # Write PETGEM connectivity in PETSc format
            writeParallelDenseMatrix(out_path, matrix, communicator=PETSc.COMM_SELF)
            # Remove temporal matrix
            del matrix

            # ---------------------------------------------------------------
            # Preprocessing edges connectivity
            # ---------------------------------------------------------------
            Print.master('     Edges connectivity')

            # Compute edges
            elemsE, edgesNodes = computeEdges(mesh.cells[0][1][:], nElems)
            nEdges = edgesNodes.shape[0]

            # Get matrix dimensions
            size = elemsE.shape

            # Build PETSc structures
            matrix = createSequentialDenseMatrixWithArray(size[0], size[1], elemsE)

            # Build path to save the file
            out_path = out_dir + '/edges.dat'

            # Write PETGEM edges in PETSc format
            writeParallelDenseMatrix(out_path, matrix, communicator=PETSc.COMM_SELF)
            # Remove temporal matrix
            del matrix

            # Reshape edgesNodes and save
            num_nodes_per_edge = 2
            num_edges_per_element = 6
            data = np.array((edgesNodes[elemsE[:], :]), dtype=np.float)
            data = data.reshape(nElems, num_nodes_per_edge*num_edges_per_element)

            # Get matrix dimensions
            size = data.shape

            # Build PETSc structures
            matrix = createSequentialDenseMatrixWithArray(size[0], size[1], data)

            # Build path to save the file
            out_path = out_dir + '/edgesNodes.dat'

            # Write PETGEM edgesNodes in PETSc format
            writeParallelDenseMatrix(out_path, matrix, communicator=PETSc.COMM_SELF)
            # Remove temporal matrix
            del matrix

            # ---------------------------------------------------------------
            # Preprocessing faces connectivity
            # ---------------------------------------------------------------
            Print.master('     Faces connectivity')

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

            # Get matrix dimensions
            size = elemsF.shape

            # Build PETSc structures
            matrix = createSequentialDenseMatrixWithArray(size[0], size[1], elemsF)

            # Build path to save the file
            out_path = out_dir + '/faces.dat'

            # Write PETGEM edges in PETSc format
            writeParallelDenseMatrix(out_path, matrix, communicator=PETSc.COMM_SELF)
            # Remove temporal matrix
            del matrix

            # ---------------------------------------------------------------
            # Preprocessing faces-edges connectivity
            # ---------------------------------------------------------------
            Print.master('     Faces-edges connectivity')

            facesE = computeFacesEdges(elemsF, elemsE, nFaces, nElems)

            num_faces_per_element = 4
            num_edges_per_face = 3
            data = np.array((facesE[elemsF[:], :]), dtype=np.float)
            data = data.reshape(nElems, num_faces_per_element*num_edges_per_face)

            # Get matrix dimensions
            size = data.shape

            # Build PETSc structures
            matrix = createSequentialDenseMatrixWithArray(size[0], size[1], data)

            # Build path to save the file
            out_path = out_dir + '/facesEdges.dat'

            # Write PETGEM edges in PETSc format
            writeParallelDenseMatrix(out_path, matrix, communicator=PETSc.COMM_SELF)
            del matrix

            # ---------------------------------------------------------------
            # Preprocessing dofs connectivity
            # ---------------------------------------------------------------
            Print.master('     DOFs connectivity')

            # Compute degrees of freedom connectivity
            basis_order = run.get('nord')
            dofs, dof_edges, dof_faces, _, total_num_dofs = computeConnectivityDOFS(elemsE,elemsF,basis_order)

            # Get matrix dimensions
            size = dofs.shape

            # Build PETSc structures
            matrix = createSequentialDenseMatrixWithArray(size[0], size[1], dofs)

            # Build path to save the file
            out_path = out_dir + '/dofs.dat'

            # Write PETGEM edges in PETSc format
            writeParallelDenseMatrix(out_path, matrix, communicator=PETSc.COMM_SELF)
            del matrix

            # ---------------------------------------------------------------
            # Preprocessing sigma model
            # ---------------------------------------------------------------
            Print.master('     Conductivity model')

            i_model = data_model.get('sigma')

            if (run.get('conductivity_from_file')):
                Print.master('     Interpolation from file not supported.')
                Print.master('     Using a constant conductivity model.')
                #exit(-1)
                # Add function to interpolate data from file
                # Allocate conductivity array
                conductivityModel = np.ones((nElems, 2), dtype=np.float)
            else:
                # Get physical groups
                elemsS = mesh.cell_data['gmsh:physical'][0]
                elemsS -= np.int(1)     # 0-based indexing

                # Get horizontal sigma
                horizontal_sigma = i_model.get('horizontal')
                vertical_sigma = i_model.get('vertical')

                # Allocate conductivity array
                conductivityModel = np.zeros((nElems, 2), dtype=np.float)

                for i in np.arange(nElems):
                    # Set horizontal sigma
                    conductivityModel[i, 0] = horizontal_sigma[np.int(elemsS[i])]

                    # Set vertical sigma
                    conductivityModel[i, 1] = vertical_sigma[np.int(elemsS[i])]

            # Get matrix dimensions
            size = conductivityModel.shape

            # Build PETSc structures
            matrix = createSequentialDenseMatrixWithArray(size[0], size[1], conductivityModel)

            # Build path to save the file
            out_path = out_dir + '/conductivityModel.dat'

            # Write PETGEM edges in PETSc format
            writeParallelDenseMatrix(out_path, matrix, communicator=PETSc.COMM_SELF)
            del matrix

            # ---------------------------------------------------------------
            # Preprocessing boundaries
            # ---------------------------------------------------------------
            Print.master('     Boundaries')

            # Compute boundary faces
            bFacesN, bFaces, nbFaces = computeBoundaryFaces(elemsF, facesN)

            # Build array with boundary dofs for csem mode (dirichlet BC)
            if (mode == 'csem'):
                # Compute boundary edges
                bEdges = computeBoundaryEdges(edgesNodes, bFacesN)

                # Compute dofs on boundaries
                _, indx_boundary_dofs = computeBoundaries(dofs, dof_edges, dof_faces, bEdges, bFaces, basis_order);

                # Build PETSc structures
                vector = createSequentialVectorWithArray(indx_boundary_dofs)

                # Build path to save the file
                out_path = out_dir + '/boundaries.dat'

                # Write PETGEM nodes in PETSc format
                writePetscVector(out_path, vector, communicator=PETSc.COMM_SELF)
                del vector

            elif (mode == 'mt'):
                # Compute to what plane the boundary face belongs
                planeFace = computeFacePlane(mesh.points, bFaces, bFacesN)

                # Compute boundary elements
                bElems, numbElems = computeBoundaryElements(elemsF, bFaces, nFaces)

                if (nbFaces != numbElems):
                    Print.master('     Number of boundary faces is not consistent.')
                    exit(-1)

                # Allocate
                data_boundaries = np.zeros((nbFaces, 53+num_dof_in_element), dtype=np.float)

                # Fill tmp matrix with data for boundary faces
                for i in np.arange(nbFaces):
                    # Get index of tetrahedral element (boundary element)
                    iEle = bElems[i]
                    # Get dofs of element container
                    dofsElement = dofs[iEle, :]

                    # Get indexes of nodes for i-boundary element and insert
                    nodesBoundaryElement = mesh.cells[0][1][iEle,:]
                    data_boundaries[i, 0:4] = nodesBoundaryElement
                    # Get nodes coordinates for i-boundary element and insert
                    coordEle = mesh.points[nodesBoundaryElement, :]
                    coordEle = coordEle.flatten()
                    data_boundaries[i, 4:16] = coordEle
                    # Get indexes of faces for i-boundary element and insert
                    facesBoundaryElement = elemsF[iEle, :]
                    data_boundaries[i, 16:20] = facesBoundaryElement
                    # Get edges indexes for faces in i-boundary element and insert
                    edgesBoundaryFace = facesE[facesBoundaryElement, :]
                    edgesBoundaryFace = edgesBoundaryFace.flatten()
                    data_boundaries[i, 20:32] = edgesBoundaryFace
                    # Get indexes of edges for i-boundary and insert
                    edgesBoundaryElement = elemsE[iEle, :]
                    data_boundaries[i, 32:38] = edgesBoundaryElement
                    # Get node indexes for edges in i-boundary and insert
                    edgesNodesBoundaryElement = edgesNodes[edgesBoundaryElement, :]
                    edgesNodesBoundaryElement = edgesNodesBoundaryElement.flatten()
                    data_boundaries[i, 38:50] = edgesNodesBoundaryElement
                    # Get plane face
                    ifacetype = planeFace[i]
                    data_boundaries[i, 50] = ifacetype
                    # Get global face index
                    localFaceIndex = bFaces[i]
                    data_boundaries[i, 51] = localFaceIndex
                    # Get sigma value
                    sigmaEle = conductivityModel[iEle, 0]
                    data_boundaries[i, 52] = sigmaEle
                    # Get dofs for boundary element and insert
                    dofsBoundaryElement = dofsElement
                    data_boundaries[i, 53::] = dofsBoundaryElement

                # Get matrix dimensions
                size = data_boundaries.shape

                # Build PETSc structures
                matrix = createSequentialDenseMatrixWithArray(size[0], size[1], data_boundaries)

                # Build path to save the file
                out_path = out_dir + '/boundaryElements.dat'

                # Write PETGEM receivers in PETSc format
                writeParallelDenseMatrix(out_path, matrix, communicator=PETSc.COMM_SELF)
                del matrix
                del data_boundaries


            # ---------------------------------------------------------------
            # Preprocessing receivers
            # ---------------------------------------------------------------
            Print.master('     Receivers')

            # 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]

            # Find out which tetrahedral element source point is in (only for csem mode)
            if (mode == 'csem'):
                # Allocate vector to save source data
                data_source = np.zeros(50+num_dof_in_element, dtype=np.float)

                i_model = data_model.get('source')

                # Get source position
                i_source_position = np.asarray(i_model.get('position'), dtype=np.float)

                # Build Delaunay triangulation with nodes
                tri = Delaunay(mesh.points)

                # Overwrite Delaunay structure with mesh_file connectivity and points
                tri.simplices = mesh.cells[0][1][:].astype(np.int32)
                tri.vertices = mesh.cells[0][1][:].astype(np.int32)

                srcElem = tri.find_simplex(i_source_position, bruteforce=True, tol=1.e-12)

                # If srcElem=-1, source not located
                if srcElem < 0:
                    Print.master('        Source no located in the computational domain. Please, verify source position or improve the mesh quality.')
                    exit(-1)

                # Build data for source insertion
                # Get indexes of nodes for srcElem and insert
                nodesSource = mesh.cells[0][1][srcElem,:]
                data_source[0:4] = nodesSource
                # Get nodes coordinates for srcElem and insert
                coordSource = mesh.points[nodesSource, :]
                coordSource = coordSource.flatten()
                data_source[4:16] = coordSource
                # Get indexes of faces for srcElem and insert
                facesSource = elemsF[srcElem, :]
                data_source[16:20] = facesSource
                # Get edges indexes for faces in srcElem and insert
                edgesFace = facesE[facesSource, :]
                edgesFace = edgesFace.flatten()
                data_source[20:32] = edgesFace
                # Get indexes of edges for srcElem and insert
                edgesSource = elemsE[srcElem, :]
                data_source[32:38] = edgesSource
                # Get node indexes for edges in srcElem and insert
                edgesNodesSource = edgesNodes[edgesSource, :]
                edgesNodesSource = edgesNodesSource.flatten()
                data_source[38:50] = edgesNodesSource
                # Get dofs for srcElem and insert
                dofsSource = dofs[srcElem,:]
                data_source[50::] = dofsSource

                # Get matrix dimensions
                size = data_source.shape

                # Build PETSc structures
                vector = createSequentialVectorWithArray(data_source)

                # Build path to save the file
                out_path = out_dir + '/source.dat'

                # Write PETGEM nodes in PETSc format
                writePetscVector(out_path, vector, communicator=PETSc.COMM_SELF)
                del vector

            # ---------------------------------------------------------------
            # Sparsity pattern
            # ---------------------------------------------------------------
            # Setup valence for each basis order (adding a small percentage to keep safe)
            valence = np.array([50, 200, 400, 800, 1400, 2500])

            # Build nnz pattern for each row
            nnz = np.full((total_num_dofs), valence[basis_order-1], dtype=np.int)

            # Build PETSc structures
            vector = createSequentialVectorWithArray(nnz)

            # Build path to save the file
            out_path = out_dir + '/nnz.dat'

            # Write PETGEM nodes in PETSc format
            writePetscVector(out_path, vector, communicator=PETSc.COMM_SELF)

            # ---------------------------------------------------------------
            # Print mesh statistics
            # ---------------------------------------------------------------
            Print.master(' ')
            Print.master('  Mesh statistics')
            Print.master('     Mesh file:            {0:12}'.format(str(model.get('mesh'))))
            Print.master('     Number of elements:   {0:12}'.format(str(nElems)))
            Print.master('     Number of faces:      {0:12}'.format(str(nFaces)))
            Print.master('     Number of edges:      {0:12}'.format(str(nEdges)))
            Print.master('     Number of dofs:       {0:12}'.format(str(total_num_dofs)))
            if (mode == 'csem'):
                Print.master('     Number of boundaries: {0:12}'.format(str(len(indx_boundary_dofs))))

            # ---------------------------------------------------------------
            # Print data model
            # ---------------------------------------------------------------
            Print.master(' ')
            Print.master('  Model data')
            Print.master('     Modeling mode:       {0:12}'.format(str(mode)))
            i_sigma = data_model.get('sigma')

            if (run.get('conductivity_from_file')):
                Print.master('     Conductivity file:   {0:12}'.format(i_sigma.get('file')))
            else:
                Print.master('     Horizontal conductivity:  {0:12}'.format(str(i_sigma.get('horizontal'))))
                Print.master('     Vertical conductivity:    {0:12}'.format(str(i_sigma.get('vertical'))))

            if (mode == 'csem'):
                i_source = data_model.get('source')
                Print.master('     Source:')
                Print.master('      - Frequency (Hz):  {0:12}'.format(str(i_source.get('frequency'))))
                Print.master('      - Position (xyz):  {0:12}'.format(str(i_source.get('position'))))
                Print.master('      - Azimuth:         {0:12}'.format(str(i_source.get('azimuth'))))
                Print.master('      - Dip:             {0:12}'.format(str(i_source.get('dip'))))
                Print.master('      - Current:         {0:12}'.format(str(i_source.get('current'))))
                Print.master('      - Length:          {0:12}'.format(str(i_source.get('length'))))
            else:
                Print.master('     Frequency (Hz):           {0:12}'.format(str(data_model.get('frequency'))))
                Print.master('     Polarization:             {0:12}'.format(str(data_model.get('polarization'))))

            Print.master('     Vector basis order:       {0:12}'.format(str(basis_order)))
            Print.master('     Receivers file:           {0:12}'.format(str(model.get('receivers'))))
            Print.master('     Number of receivers:      {0:12}'.format(str(nReceivers)))
            Print.master('     VTK output:               {0:12}'.format(str(output.get('vtk'))))
            Print.master('     Cuda support:             {0:12}'.format(str(run.get('cuda'))))
            Print.master('     Output directory:         {0:12}'.format(str(output.get('directory'))))
            Print.master('     Scratch directory:        {0:12}'.format(str(output.get('directory_scratch'))))

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

        # Apply barrier for MPI tasks alignement
        parEnv.comm.barrier()

        return

# ###############################################################
# ################     FUNCTIONS DEFINITION     #################
# ###############################################################
[docs]def unitary_test(): """Unitary test for preprocessing.py script."""
# ############################################################### # ################ MAIN ################# # ############################################################### if __name__ == '__main__': # --------------------------------------------------------------- # Run unitary test # --------------------------------------------------------------- unitary_test()