Pre/Post deep learning

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

alternate text

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):

  1import numpy as np
  2
  3#################
  4# PREPROCESSING #
  5#################
  6
  7# Read the solution generated by a physical code
  8from Muscat.IO import XdmfReader as XR
  9reader = XR.XdmfReader(filename = "PrePostDeepLearning_Input.xmf")
 10reader.Read()
 11
 12dom = reader.xdmf.GetDomain(0)
 13grid = dom.GetGrid(0)
 14
 15# Read the mesh
 16uMesh = grid.GetSupport()
 17
 18# Read the solution field 'U'
 19indexU = grid.GetPointFieldsNames().index("U")
 20U = grid.GetPointFields()[indexU][:,0:2]
 21
 22# Create a rectilinear mesh of size 48*48 excluding the left part of the mesh (x<0)
 23Nx = 48
 24Ny = 48
 25boundingMin, boundingMax = uMesh.ComputeBoundingBox()
 26Lx = boundingMax[0] - 0.
 27Ly = boundingMax[1] - boundingMin[1]
 28
 29from Muscat.Containers.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh
 30unstructuredRectMesh = CreateConstantRectilinearMesh(dimensions=[Nx,Ny], spacing=[Lx/(Nx-1), Ly/(Ny-1)], origin=[0., boundingMin[1]])
 31
 32# Compute the projection operator from unstructured mesh to structured mesh
 33from Muscat.FE.FETools import PrepareFEComputation
 34from Muscat.Containers.MeshFieldOperations import GetFieldTransferOp
 35from Muscat.FE.Fields.FEField import FEField
 36
 37space, numberings, _offset, _NGauss = PrepareFEComputation(uMesh, numberOfComponents = 2)
 38inputFEField = FEField(name="U",mesh=uMesh, space=space, numbering=numberings[0])
 39methods = ["Interp/Nearest","Nearest/Nearest","Interp/Clamp","Interp/Extrap","Interp/ZeroFill"]
 40operator, status = GetFieldTransferOp(inputFEField, unstructuredRectMesh.nodes, method = methods[2], verbose=True)
 41
 42# Compute the projected field on the structured mesh
 43projectedU = operator.dot(U)
 44
 45# Export the structured mesh and projected field in xdmf format (in ASCII)
 46# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)
 47from Muscat.IO import XdmfWriter as XW
 48writer = XW.XdmfWriter('PrePostDeepLearning_OutputI.xdmf')
 49writer.SetHdf5(False)
 50writer.Open()
 51writer.Write(unstructuredRectMesh,PointFields=[projectedU], PointFieldsNames=["U"])
 52writer.Close()
 53
 54
 55##################
 56# POSTPROCESSING #
 57##################
 58
 59# Compute the projection operator from the structured mesh to the unstructured mesh (inverse projection)
 60space, numberings, _offset, _NGauss = PrepareFEComputation(unstructuredRectMesh)
 61inputFEField = FEField(name="U",mesh=unstructuredRectMesh,space=space,numbering=numberings[0])
 62operator, status = GetFieldTransferOp(inputFEField, uMesh.nodes, method = methods[2], verbose=True)
 63
 64# Compute the inverse-projected projected field on the unstructured mesh
 65inverseProjected_ProjectedU = operator.dot(projectedU)
 66
 67# Export the unstructured mesh and inverse-projected projected field in xdmf format
 68# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)
 69writer = XW.XdmfWriter('PrePostDeepLearning_OutputII.xdmf')
 70writer.SetHdf5(False)
 71writer.Open()
 72writer.Write(uMesh, PointFields=[inverseProjected_ProjectedU], PointFieldsNames=["U"])
 73writer.Close()
 74
 75
 76##################
 77# CHECK ACCURACY #
 78##################
 79
 80# Create a clippedMesh from the unstructured mesh by removing the left part (x<0)
 81# (with field, inverse-projected projected field and difference between them)
 82from Muscat.Containers.Filters.FilterObjects import ElementFilter
 83from Muscat.FE.FETools import ComputePhiAtIntegPoint, ComputeL2ScalarProducMatrix
 84from Muscat.Containers.MeshInspectionTools import ExtractElementsByElementFilter
 85from Muscat.Containers.MeshModificationTools import CleanLonelyNodes
 86from Muscat.Containers.MeshFieldOperations import CopyFieldsFromOriginalMeshToTargetMesh
 87
 88uMesh.nodeFields['delta'] = U - inverseProjected_ProjectedU
 89
 90elFilter = ElementFilter(zone = lambda p: (-p[:,0]))
 91meshClipped = ExtractElementsByElementFilter(uMesh, elFilter)
 92CleanLonelyNodes(meshClipped)
 93
 94CopyFieldsFromOriginalMeshToTargetMesh(uMesh, meshClipped)
 95nbeOfNodes = meshClipped.GetNumberOfNodes()
 96
 97deltaClippedMesh = meshClipped.nodeFields['delta']
 98
 99from Muscat.Helpers.TextFormatHelper import TFormat
100from Muscat.Helpers.Timer import Timer
101
102print("Compute the L2(Omega) norm of delta by applying numerical quadrature from Lagrange P1")
103print("Finite element integration using three different methods")
104
105#### Method 1
106print(TFormat.Center(TFormat.InRed("Method 1: ")+TFormat.InBlue("'by hand'")))
107
108timer = Timer("Duration of method 1")
109#compute method 1 three times
110for i in range(3):
111    timer.Start()
112    integrationWeights, phiAtIntegPoint = ComputePhiAtIntegPoint(meshClipped)
113
114    vectDeltaAtIntegPoints = np.empty((2,phiAtIntegPoint.shape[0]))
115    for i in range(2):
116        vectDeltaAtIntegPoints[i] = phiAtIntegPoint.dot(deltaClippedMesh[:,i])
117
118    normDeltaMethod1 = np.sqrt(np.einsum('ij,ij,j', vectDeltaAtIntegPoints, vectDeltaAtIntegPoints, integrationWeights, optimize = True))
119    timer.Stop()
120
121print("norm(Delta) =", normDeltaMethod1)
122############
123
124#### Method 2
125print(TFormat.Center(TFormat.InRed("Method 2: ")+TFormat.InBlue("using the mass matrix")))
126
127timer = Timer("Duration of method 2").Start()
128
129massMatrix = ComputeL2ScalarProducMatrix(meshClipped, 2)
130
131normDeltaMethod2 = np.sqrt(np.dot(deltaClippedMesh.ravel(order='F'), massMatrix.dot(deltaClippedMesh.ravel(order='F'))))
132
133timer.Stop()
134print("norm(Delta) =", normDeltaMethod2)
135############
136
137#### Method 3
138print(TFormat.Center(TFormat.InRed("Method 3: ")+TFormat.InBlue("using the weak form engine")))
139
140from Muscat.FE.Integration import IntegrateGeneral
141from Muscat.FE.SymWeakForm import GetField, GetTestField
142from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1, ConstantSpaceGlobal
143from Muscat.FE.DofNumbering import ComputeDofNumbering
144
145timer = Timer("Duration of method 3").Start()
146
147U = GetField("U",2)
148Tt = GetTestField("T",1)
149
150wf = U.T*U*Tt
151
152#meshClipped
153numbering = ComputeDofNumbering(meshClipped,LagrangeSpaceP1)
154field1 = FEField("U_0",mesh=meshClipped,space=LagrangeSpaceP1,numbering=numbering, data = deltaClippedMesh[:,0])
155field1.ConvertDataForNativeTreatment()
156field2 = FEField("U_1",mesh=meshClipped,space=LagrangeSpaceP1,numbering=numbering, data = deltaClippedMesh[:,1])
157field2.ConvertDataForNativeTreatment()
158
159numbering = ComputeDofNumbering(meshClipped,ConstantSpaceGlobal)
160unknownField = FEField("T",mesh=meshClipped,space=ConstantSpaceGlobal,numbering=numbering)
161
162elFilter = ElementFilter()
163
164K, F = IntegrateGeneral(mesh=meshClipped,
165                    wform=wf,
166                    constants={},
167                    fields=[field1,field2],
168                    unknownFields=[unknownField],
169                    integrationRuleName="LagrangeP1Quadrature",
170                    elementFilter=elFilter)
171
172normDeltaMethod3 = np.sqrt(F[0])
173
174timer.Stop()
175print("norm(Delta) =", normDeltaMethod2)
176############
177
178print(Timer.PrintTimes())
179
180
181