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.

alternate text

Fig. 1 Example of deep learning pre/post

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

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