{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a6f88c8a-c1cf-4923-94fd-99362247caa4",
   "metadata": {},
   "source": [
    "# Pre/Post deep learning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e17c05df-a6e7-4ba1-8475-323718ca3e4e",
   "metadata": {},
   "source": [
    "Some deep learning workflows applied to physics contexts require the projection of fields defined on an unstructured mesh onto a rectilinear grid, and inversely."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c90b07d3-ae47-4a78-8e51-2404c82d518e",
   "metadata": {},
   "source": [
    "## Includes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "065452c6-4d36-493d-ac2a-6cf10f5f49a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyvista\n",
    "pyvista.global_theme._jupyter_backend = 'panel' # remove this line to get interactive 3D plots\n",
    "import numpy as np\n",
    "from Muscat.Bridges.PyVistaBridge import PlotMesh, MeshToPyVista\n",
    "from Muscat.Helpers.IO.TemporaryDirectory import TemporaryDirectory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "129cd9af-74c1-44c2-96c1-883d2aca4fd0",
   "metadata": {},
   "source": [
    "## PREPROCESSING"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f87c4d8b-65d3-4ec0-bbe7-f3e250feeba3",
   "metadata": {},
   "source": [
    "### Read the solution generated by a physical code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd3da3c5-67f8-4ff9-ac52-c4125ab90b77",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Muscat.IO import XdmfReader as XR\n",
    "reader = XR.XdmfReader(filename = \"PrePostDeepLearning_Input.xmf\")\n",
    "reader.Read()\n",
    "grid = reader.xdmf.GetDomain(0).GetGrid(0)\n",
    "\n",
    "# Read the mesh\n",
    "uMesh = grid.GetSupport()\n",
    "\n",
    "# Read the solution field 'U'\n",
    "indexU = grid.GetPointFieldsNames().index(\"U\")\n",
    "U = grid.GetPointFields()[indexU][:,0:2]\n",
    "uMesh.nodeFields[\"U\"] = U\n",
    "PlotMesh(uMesh,scalars=\"U\", show_edges=True, cpos=\"xy\", show_scalar_bar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96c300ff-575d-4832-bfe4-f07a0e0425dd",
   "metadata": {},
   "source": [
    "### Create a rectilinear mesh of size 48*48 excluding the left part of the mesh (x<0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "170d1c64-9e68-4e7d-b9b9-91ef441a5027",
   "metadata": {},
   "outputs": [],
   "source": [
    "Nx = 48\n",
    "Ny = 48\n",
    "boundingMin, boundingMax = uMesh.ComputeBoundingBox()\n",
    "Lx = boundingMax[0] - 0.\n",
    "Ly = boundingMax[1] - boundingMin[1]\n",
    "\n",
    "from Muscat.MeshTools.ConstantRectilinearMeshTools import CreateConstantRectilinearMesh\n",
    "unstructuredRectMesh = CreateConstantRectilinearMesh(dimensions=[Nx,Ny], spacing=[Lx/(Nx-1), Ly/(Ny-1)], origin=[0., boundingMin[1]])\n",
    "PlotMesh(unstructuredRectMesh, show_edges=True, cpos=\"xy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59e3e1e9-5c07-4f69-a21d-817cb84a21de",
   "metadata": {},
   "source": [
    "## Compute the projection operator from unstructured mesh to structured mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c54fcbb-88ac-43b2-b0d4-aa6ddc81c629",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Muscat.FE.FETools import PrepareFEComputation\n",
    "from Muscat.MeshTools.MeshFieldOperations import GetFieldTransferOp\n",
    "from Muscat.FE.Fields.FEField import FEField\n",
    "\n",
    "inputFEField = FEField(name=\"U\",mesh=uMesh)\n",
    "\n",
    "operator, status, entities  = GetFieldTransferOp(inputFEField, unstructuredRectMesh.nodes, method = \"Interp/Clamp\", verbose=True)\n",
    "\n",
    "# Compute the projected field on the structured mesh\n",
    "projectedU = operator.dot(U)\n",
    "\n",
    "unstructuredRectMesh.nodeFields[\"projectedU\"] = projectedU\n",
    "PlotMesh(unstructuredRectMesh,scalars=\"projectedU\", show_edges=True, cpos=\"xy\", show_scalar_bar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75188b5f-b0bc-4217-86b5-910acb6c623e",
   "metadata": {},
   "source": [
    "## Export the structured mesh and projected field in xdmf format (in ASCII)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1aaf926e-23f9-4ad3-a8b1-22c1702577b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To visualize this xdmf file you can use ParaView (downloadable from https://www.paraview.org/)\n",
    "from Muscat.IO import XdmfWriter as XW\n",
    "writer = XW.XdmfWriter(TemporaryDirectory.GetTempPath() + 'PrePostDeepLearning_OutputI.xdmf')\n",
    "writer.SetHdf5(False)\n",
    "writer.Open()\n",
    "writer.Write(unstructuredRectMesh,PointFields=[projectedU], PointFieldsNames=[\"U\"])\n",
    "writer.Close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc4837a0-105c-465b-8c33-cbfb1c4e6375",
   "metadata": {},
   "source": [
    "## POSTPROCESSING"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e94f9c3b-cae6-4caa-87e0-1c9e86054817",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the projection operator from the structured mesh to the unstructured mesh (inverse projection)\n",
    "inputFEField = FEField(name=\"U\",mesh=unstructuredRectMesh)\n",
    "operator, status, entities = GetFieldTransferOp(inputFEField, uMesh.nodes, method =  \"Interp/Clamp\", verbose=True)\n",
    "\n",
    "# Compute the inverse-projected projected field on the unstructured mesh\n",
    "inverseProjected_ProjectedU = operator.dot(projectedU)\n",
    "\n",
    "uMesh.nodeFields[\"inverseProjected_ProjectedU\"] = inverseProjected_ProjectedU\n",
    "PlotMesh(uMesh,scalars=\"inverseProjected_ProjectedU\", show_edges=False, cpos=\"xy\", show_scalar_bar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "248e4af1-a4dd-4c3a-888a-3ed3873b251c",
   "metadata": {},
   "source": [
    "## CHECK ACCURACY"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56a8b833-1488-4d31-b6f6-5c0d56324bdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a clippedMesh from the unstructured mesh by removing the left part (x<0)\n",
    "# (with field, inverse-projected projected field and difference between them)\n",
    "from Muscat.MeshContainers.Filters.FilterObjects import ElementFilter\n",
    "from Muscat.FE.FETools import ComputePhiAtIntegPoint, ComputeL2ScalarProductMatrix\n",
    "from Muscat.MeshTools.MeshInspectionTools import ExtractElementsByElementFilter\n",
    "from Muscat.MeshTools.MeshModificationTools import CleanLonelyNodes\n",
    "from Muscat.MeshTools.MeshFieldOperations import CopyFieldsFromOriginalMeshToTargetMesh\n",
    "\n",
    "uMesh.nodeFields['delta'] = U - inverseProjected_ProjectedU\n",
    "\n",
    "elFilter = ElementFilter(zone = lambda p: (-p[:,0]))\n",
    "meshClipped = ExtractElementsByElementFilter(uMesh, elFilter)\n",
    "CleanLonelyNodes(meshClipped)\n",
    "\n",
    "CopyFieldsFromOriginalMeshToTargetMesh(uMesh, meshClipped)\n",
    "nbeOfNodes = meshClipped.GetNumberOfNodes()\n",
    "\n",
    "deltaClippedMesh = meshClipped.nodeFields['delta']\n",
    "\n",
    "PlotMesh(meshClipped,scalars=\"delta\", show_edges=False, cpos=\"xy\", show_scalar_bar=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "257b31f3-f60f-42b1-bc46-360ce2259144",
   "metadata": {},
   "source": [
    "## Computation of the integral of the delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56947362-9398-4ecd-aeb5-4112296a7867",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Muscat.Helpers.TextFormatHelper import TFormat\n",
    "from Muscat.Helpers.Timer import Timer\n",
    "\n",
    "print(\"Compute the L2(Omega) norm of delta by applying numerical quadrature from Lagrange P1\")\n",
    "print(\"Finite element integration using three different methods\")\n",
    "\n",
    "#### Method 1 #######################################################################################################################\n",
    "print(TFormat.Center(TFormat.InRed(\"Method 1: \")+TFormat.InBlue(\"'by hand'\")))\n",
    "\n",
    "timer = Timer(\"Duration of method 1\")\n",
    "#compute method 1 three times\n",
    "for i in range(2):\n",
    "    timer.Start()\n",
    "    integrationWeights, phiAtIntegPoint = ComputePhiAtIntegPoint(meshClipped)\n",
    "\n",
    "    vectDeltaAtIntegPoints = np.empty((2,phiAtIntegPoint.shape[0]))\n",
    "    for i in range(2):\n",
    "        vectDeltaAtIntegPoints[i] = phiAtIntegPoint.dot(deltaClippedMesh[:,i])\n",
    "\n",
    "    normDeltaMethod1 = np.sqrt(np.einsum('ij,ij,j', vectDeltaAtIntegPoints, vectDeltaAtIntegPoints, integrationWeights, optimize = True))\n",
    "    timer.Stop()\n",
    "\n",
    "print(\"norm(Delta) =\", normDeltaMethod1)\n",
    "############\n",
    "#### Method 2 #######################################################################################################################\n",
    "print(TFormat.Center(TFormat.InRed(\"Method 2: \")+TFormat.InBlue(\"using the mass matrix\")))\n",
    "\n",
    "timer = Timer(\"Duration of method 2\").Start()\n",
    "\n",
    "massMatrix = ComputeL2ScalarProductMatrix(meshClipped, 2)\n",
    "\n",
    "normDeltaMethod2 = np.sqrt(np.dot(deltaClippedMesh.ravel(order='F'), massMatrix.dot(deltaClippedMesh.ravel(order='F'))))\n",
    "\n",
    "timer.Stop()\n",
    "print(\"norm(Delta) =\", normDeltaMethod2)\n",
    "\n",
    "#### Method 3 #######################################################################################################################\n",
    "print(TFormat.Center(TFormat.InRed(\"Method 3: \")+TFormat.InBlue(\"using Muscat integration tool\")))\n",
    "print(TFormat.Center(TFormat.InRed(\"Warning \")+TFormat.InBlue(\"this method squares the values at the nodes before integrating: not recommended\")))\n",
    "\n",
    "from Muscat.FE.IntegrationTools import IntegrateField\n",
    "timer = Timer(\"Duration of method 3\").Start()\n",
    "deltaSquare_U0 = FEField(\"delta_U0\",meshClipped,data= deltaClippedMesh[:,0].flatten() )**2 + FEField(\"delta_U1\",meshClipped,data= deltaClippedMesh[:,1].flatten() )**2\n",
    "F = IntegrateField(deltaSquare_U0, ElementFilter(dimensionality=2))\n",
    "timer.Stop()\n",
    "print(\"norm(Delta) =\", np.sqrt(F))\n",
    "\n",
    "#### Method 4 #######################################################################################################################\n",
    "print(TFormat.Center(TFormat.InRed(\"Method 4: \")+TFormat.InBlue(\"using the weak form engine\")))\n",
    "\n",
    "from Muscat.FE.Integration import IntegrateGeneral\n",
    "from Muscat.FE.SymWeakForm import GetField, GetTestField\n",
    "from Muscat.FE.Spaces.FESpaces import LagrangeSpaceP1, ConstantSpaceGlobal\n",
    "from Muscat.FE.DofNumbering import ComputeDofNumbering\n",
    "\n",
    "timer = Timer(\"Duration of method 4\").Start()\n",
    "\n",
    "U = GetField(\"U\",2,2)\n",
    "Tt = GetTestField(\"T\",1,2)\n",
    "\n",
    "wf = U.T*U*Tt\n",
    "\n",
    "#meshClipped\n",
    "\n",
    "field1 = FEField(\"U_0\",mesh=meshClipped, data = deltaClippedMesh[:,0].flatten())\n",
    "field2 = FEField(\"U_1\",mesh=meshClipped, data = deltaClippedMesh[:,1].flatten())\n",
    "\n",
    "numbering = ComputeDofNumbering(meshClipped,ConstantSpaceGlobal)\n",
    "unknownField = FEField(\"T\",mesh=meshClipped,space=ConstantSpaceGlobal,numbering=numbering)\n",
    "\n",
    "elFilter = ElementFilter()\n",
    "\n",
    "K, F = IntegrateGeneral(mesh=meshClipped,\n",
    "                    wform=wf,\n",
    "                    constants={},\n",
    "                    fields=[field1,field2],\n",
    "                    unknownFields=[unknownField],\n",
    "                    integrationRuleName=\"LagrangeP1Quadrature\",\n",
    "                    elementFilter=elFilter)\n",
    "\n",
    "normDeltaMethod3 = np.sqrt(F[0])\n",
    "\n",
    "timer.Stop()\n",
    "print(\"norm(Delta) =\", normDeltaMethod2)\n",
    "############\n",
    "\n",
    "print(Timer.PrintTimes())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}