Pre/Post deep learning

Some deep learning workflow applied to physics contexts require the projection of fields defined on an unstructured mesh onto a rectilinear grid, and inversely.

Fig. 2 Example of deep learning pre/post

The mesh and solution associated to a previously computed flow field is read (see also Fig. 2, top-left image):

import numpy as np

#################
# PREPROCESSING #
#################

# Read the solution generated by a physical code
from BasicTools.IO import XdmfReader as XR
reader = XR.XdmfReader(filename = "PrePostDeepLearning_Input.xmf")
reader.Read()

dom = reader.xdmf.GetDomain(0)
grid = dom.GetGrid(0)

# Read the mesh
uMesh = grid.GetSupport()
# convert all data to the correct binary (float64, int64) representation for the c++ part of BasicTools
uMesh.ConvertDataForNativeTreatment()

# Read the solution field 'U'
indexU = grid.GetPointFieldsNames().index("U")
U = grid.GetPointFields()[indexU][:,0:2]

# Create a rectilinear mesh of size 48*48 excluding the left part of the mesh (x<0)
Nx = 48
Ny = 48
uMesh.ComputeBoundingBox()
Lx = uMesh.boundingMax[0] - 0.
Ly = uMesh.boundingMax[1] - uMesh.boundingMin[1]

from BasicTools.Containers.ConstantRectilinearMesh import ConstantRectilinearMesh
rectMesh = ConstantRectilinearMesh(dim=2)
rectMesh.SetDimensions([Nx,Ny])
rectMesh.SetSpacing([Lx/(Nx-1), Ly/(Ny-1)])
rectMesh.SetOrigin([0., uMesh.boundingMin[1]])

# Compute the projection operator from unstructured mesh to structured mesh
from BasicTools.Containers.UnstructuredMeshCreationTools import CreateMeshFromConstantRectilinearMesh
from BasicTools.FE.FETools import PrepareFEComputation
from BasicTools.Containers.UnstructuredMeshFieldOperations import GetFieldTransferOp
from BasicTools.FE.Fields.FEField import FEField

unstructuredRectMesh = CreateMeshFromConstantRectilinearMesh(rectMesh)
space, numberings, _offset, _NGauss = PrepareFEComputation(uMesh, numberOfComponents = 2)
inputFEField = FEField(name="U",mesh=uMesh, space=space, numbering=numberings[0])
methods = ["Interp/Nearest","Nearest/Nearest","Interp/Clamp","Interp/Extrap","Interp/ZeroFill"]
operator, status = GetFieldTransferOp(inputFEField, unstructuredRectMesh.nodes, method = methods[2], verbose=True)

# Compute the projected field on the structured mesh
projectedU = operator.dot(U)

# Export the structured mesh and projected field in xdmf format (in ASCII)
# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)
from BasicTools.IO import XdmfWriter as XW
writer = XW.XdmfWriter('PrePostDeepLearning_OutputI.xdmf')
writer.SetHdf5(False)
writer.Open()
writer.Write(rectMesh,PointFields=[projectedU.reshape((Nx,Ny,2))], PointFieldsNames=["U"])
writer.Close()


##################
# POSTPROCESSING #
##################

# Compute the projection operator from the structured mesh to the unstructured mesh (inverse projection)
space, numberings, _offset, _NGauss = PrepareFEComputation(unstructuredRectMesh)
inputFEField = FEField(name="U",mesh=unstructuredRectMesh,space=space,numbering=numberings[0])
operator, status = GetFieldTransferOp(inputFEField, uMesh.nodes, method = methods[2], verbose=True)

# Compute the inverse-projected projected field on the unstructured mesh
inverseProjected_ProjectedU = operator.dot(projectedU)

# Export the unstructured mesh and inverse-projected projected field in xdmf format
# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)
writer = XW.XdmfWriter('PrePostDeepLearning_OutputII.xdmf')
writer.SetHdf5(False)
writer.Open()
writer.Write(uMesh, PointFields=[inverseProjected_ProjectedU], PointFieldsNames=["U"])
writer.Close()


##################
# CHECK ACCURACY #
##################

# Create a clippedMesh from the unstructured mesh by removing the left part (x<0)
# (with field, inverse-projected projected field and difference between them)
from BasicTools.Containers.Filters import ElementFilter
from BasicTools.FE.FETools import ComputePhiAtIntegPoint, ComputeL2ScalarProducMatrix
from BasicTools.Containers.UnstructuredMeshInspectionTools import ExtractElementsByElementFilter
from BasicTools.Containers.UnstructuredMeshModificationTools import CleanLonelyNodes
from BasicTools.Containers.UnstructuredMeshFieldOperations import CopyFieldsFromOriginalMeshToTargetMesh

uMesh.nodeFields['delta'] = U - inverseProjected_ProjectedU

elFilter = ElementFilter(uMesh, zone = lambda p: (-p[:,0]))
unstructuredMeshClipped = ExtractElementsByElementFilter(uMesh, elFilter)
CleanLonelyNodes(unstructuredMeshClipped)
unstructuredMeshClipped.ConvertDataForNativeTreatment()

CopyFieldsFromOriginalMeshToTargetMesh(uMesh, unstructuredMeshClipped)
nbeOfNodes = unstructuredMeshClipped.GetNumberOfNodes()

deltaClippedMesh = unstructuredMeshClipped.nodeFields['delta']

from BasicTools.Helpers.TextFormatHelper import TFormat
from BasicTools.Helpers.Timer import Timer

print("Compute the L2(Omega) norm of delta by applying numerical quadrature from Lagrange P1")
print("Finite element integration using three different methods")

#### Method 1
print(TFormat.Center(TFormat.InRed("Method 1: ")+TFormat.InBlue("'by hand'")))

timer = Timer("Duration of method 1")
#compute method 1 three times
for i in range(3):
    timer.Start()
    integrationWeights, phiAtIntegPoint = ComputePhiAtIntegPoint(unstructuredMeshClipped)

    vectDeltaAtIntegPoints = np.empty((2,phiAtIntegPoint.shape[0]))
    for i in range(2):
        vectDeltaAtIntegPoints[i] = phiAtIntegPoint.dot(deltaClippedMesh[:,i])

    normDeltaMethod1 = np.sqrt(np.einsum('ij,ij,j', vectDeltaAtIntegPoints, vectDeltaAtIntegPoints, integrationWeights, optimize = True))
    timer.Stop()

print("norm(Delta) =", normDeltaMethod1)
############

#### Method 2
print(TFormat.Center(TFormat.InRed("Method 2: ")+TFormat.InBlue("using the mass matrix")))

timer = Timer("Duration of method 2").Start()

massMatrix = ComputeL2ScalarProducMatrix(unstructuredMeshClipped, 2)

normDeltaMethod2 = np.sqrt(np.dot(deltaClippedMesh.ravel(order='F'), massMatrix.dot(deltaClippedMesh.ravel(order='F'))))

timer.Stop()
print("norm(Delta) =", normDeltaMethod2)
############

#### Method 3
print(TFormat.Center(TFormat.InRed("Method 3: ")+TFormat.InBlue("using the weak form engine")))

from BasicTools.FE.Integration import IntegrateGeneral
from BasicTools.FE.SymWeakForm import GetField, GetTestField
from BasicTools.FE.Spaces.FESpaces import LagrangeSpaceP1, ConstantSpaceGlobal
from BasicTools.FE.DofNumbering import ComputeDofNumbering

timer = Timer("Duration of method 3").Start()

U = GetField("U",2)
Tt = GetTestField("T",1)

wf = U.T*U*Tt

#unstructuredMeshClipped
numbering = ComputeDofNumbering(unstructuredMeshClipped,LagrangeSpaceP1)
field1 = FEField("U_0",mesh=unstructuredMeshClipped,space=LagrangeSpaceP1,numbering=numbering, data = deltaClippedMesh[:,0])
field1.ConvertDataForNativeTreatment()
field2 = FEField("U_1",mesh=unstructuredMeshClipped,space=LagrangeSpaceP1,numbering=numbering, data = deltaClippedMesh[:,1])
field2.ConvertDataForNativeTreatment()

numbering = ComputeDofNumbering(unstructuredMeshClipped,ConstantSpaceGlobal)
unkownField = FEField("T",mesh=unstructuredMeshClipped,space=ConstantSpaceGlobal,numbering=numbering)

elFilter = ElementFilter(unstructuredMeshClipped)

K, F = IntegrateGeneral(mesh=unstructuredMeshClipped,
                    wform=wf,
                    constants={},
                    fields=[field1,field2],
                    unkownFields=[unkownField],
                    integrationRuleName="LagrangeP1",
                    elementFilter=elFilter)

normDeltaMethod3 = np.sqrt(F[0])

timer.Stop()
print("norm(Delta) =", normDeltaMethod2)
############

print(Timer.PrintTimes())