{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is a test for the volume & surface integral functionality. \n",
"\n",
"Various mesh are integrated over, with checks made against things like total volume, surface area and function integration. "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import underworld as uw\n",
"import underworld.visualisation as vis\n",
"import numpy as np\n",
"import math"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\tGlobal element size: 16x16\n",
"\tLocal offset of rank 0: 0x0\n",
"\tLocal range of rank 0: 16x16\n",
"\tGlobal element size: 16x16\n",
"\tLocal offset of rank 0: 0x0\n",
"\tLocal range of rank 0: 16x16\n"
]
}
],
"source": [
"minCoord = [-2.,-1.]\n",
"maxCoord = [ 2., 1.]\n",
"meshq1 = uw.mesh.FeMesh_Cartesian(elementType='Q1', elementRes=(16,16), minCoord=minCoord, maxCoord=maxCoord)\n",
"meshq2 = uw.mesh.FeMesh_Cartesian(elementType='Q2', elementRes=(16,16), minCoord=minCoord, maxCoord=maxCoord)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"meshfig = vis.Figure(figsize=(800,400))\n",
"meshfig.append( vis.objects.Mesh(meshq1) )\n",
"meshfig.show()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def checkSurfaceIntegration(mesh, value, fn=1., rtol=1e-05):\n",
" \"\"\"\n",
" This function performs surface integration across the entire boundary of the \n",
" provided mesh. The surface integration of 'fn' is expected to return a result \n",
" with 'rtol' of 'value'.\n",
" \n",
" Parameters\n",
" ----------\n",
" mesh : uw.mesh.FeMesh\n",
" The mesh across which integration is to be performed.\n",
" value: float\n",
" The expected result of the integration.\n",
" fn : uw.function.Function\n",
" The function which is to be integrated.\n",
" rtol : float\n",
" Relative tolerance. Check allclose documentation for details.\n",
" \"\"\"\n",
" boundaryLength = uw.utils.Integral( fn=fn, mesh=mesh, integrationType='Surface', surfaceIndexSet=mesh.specialSets[\"AllWalls_VertexSet\"])\n",
" result = boundaryLength.evaluate()\n",
" assert np.allclose( result, value, rtol=rtol ), \"Error occurred in surface integration. Integration = {}, Expected = {}\".format(result, value)\n",
"def checkVolumeIntegration(mesh, value, fn=1., rtol=1e-05):\n",
" \"\"\"\n",
" This function performs volume integration across the entire volume of the \n",
" provided mesh. The volume integration of 'fn' is expected to return a result \n",
" with 'rtol' of 'value'.\n",
" \n",
" Parameters\n",
" ----------\n",
" mesh : uw.mesh.FeMesh\n",
" The mesh across which integration is to be performed.\n",
" value: float\n",
" The expected result of the integration.\n",
" fn : uw.function.Function\n",
" The function which is to be integrated.\n",
" rtol : float\n",
" Relative tolerance. Check allclose documentation for details.\n",
" \"\"\"\n",
"\n",
" totalVolume = uw.utils.Integral( fn=fn, mesh=mesh, integrationType='Volume')\n",
" result = totalVolume.evaluate()\n",
" assert np.allclose( result, value, rtol=rtol ), \"Error occurred in volume integration. Integration = {}, Expected = {}\".format(result, value)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def correctLength(minCoord,maxCoord):\n",
" \"\"\" Calculates the total boundary length of a mesh created using minCoord,maxCoord. \"\"\"\n",
" return 2.*sum( maxCoord[ii] - minCoord[ii] for ii in range(0,len(minCoord)))\n",
"def correctVol(minCoord,maxCoord):\n",
" \"\"\" Calculates the volume of a mesh created using minCoord,maxCoord. \"\"\"\n",
" return np.prod( [ maxCoord[ii] - minCoord[ii] for ii in range(0,len(minCoord)) ] ) "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# check we are getting correct \n",
"for mesh in [meshq1,meshq2]:\n",
" mesh.reset()\n",
" checkSurfaceIntegration(mesh, correctLength(minCoord,maxCoord)) \n",
" checkVolumeIntegration(mesh, correctVol(minCoord,maxCoord))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# lets rotate and try again\n",
"def rotate(mesh):\n",
" with mesh.deform_mesh():\n",
" for vert in mesh.data[:]:\n",
" radius = np.linalg.norm(vert[0:2])\n",
" angle = np.arctan2(vert[1],vert[0])\n",
" vert[0] = radius*math.cos(angle+math.pi/6.)\n",
" vert[1] = radius*math.sin(angle+math.pi/6.)\n",
"# check\n",
"for mesh in [meshq1,meshq2]:\n",
" mesh.reset()\n",
" rotate(mesh)\n",
" checkSurfaceIntegration(mesh, correctLength(minCoord,maxCoord)) \n",
" checkVolumeIntegration( mesh, correctVol(minCoord,maxCoord))\n",
"meshfig.show()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def refine(mesh, minCoord, maxCoord):\n",
" # lets do some refinement and check\n",
" xp = np.array( (0., 0.10, 0.20, 0.300, 0.40, 0.5, 0.60, 0.70, 0.8, 0.90, 1.0) )\n",
" yp = np.array( (0., 0.05, 0.10, 0.125, 0.15, 0.2, 0.25, 0.30, 0.4, 0.55, 1.0) )\n",
" # scale to domain\n",
" xp = xp*(maxCoord[0]-minCoord[0]) - maxCoord[0]\n",
" yp = yp*(maxCoord[0]-minCoord[0]) - maxCoord[0]\n",
" with mesh.deform_mesh():\n",
" mesh.data[:,0] = np.interp( mesh.data[:,0], xp, yp) \n",
"\n",
"# check\n",
"for mesh in [meshq1,meshq2]:\n",
" mesh.reset()\n",
" refine(mesh, minCoord, maxCoord)\n",
" checkSurfaceIntegration(mesh, correctLength(minCoord,maxCoord)) \n",
" checkVolumeIntegration( mesh, correctVol(minCoord,maxCoord))\n",
" rotate(mesh)\n",
" checkSurfaceIntegration(mesh, correctLength(minCoord,maxCoord)) \n",
" checkVolumeIntegration( mesh, correctVol(minCoord,maxCoord))\n",
"meshfig.show()"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\tGlobal element size: 16x16\n",
"\tLocal offset of rank 0: 0x0\n",
"\tLocal range of rank 0: 16x16\n",
"\tGlobal element size: 16x16\n",
"\tLocal offset of rank 0: 0x0\n",
"\tLocal range of rank 0: 16x16\n"
]
},
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# giordani mesh \n",
"minCoord = [-1.,-1.]\n",
"maxCoord = [ 1., 1.]\n",
"meshq1 = uw.mesh.FeMesh_Cartesian(elementType='Q1', elementRes=(16,16), minCoord=minCoord, maxCoord=maxCoord)\n",
"meshq2 = uw.mesh.FeMesh_Cartesian(elementType='Q2', elementRes=(16,16), minCoord=minCoord, maxCoord=maxCoord)\n",
"def circleMesh(mesh):\n",
" with mesh.deform_mesh():\n",
" for vert in mesh.data[:]:\n",
" radius = np.max(np.abs(vert))\n",
" angle = np.arctan2(vert[1],vert[0])\n",
" vert[0] = radius*math.cos(angle)\n",
" vert[1] = radius*math.sin(angle)\n",
"\n",
"# check.. note we need to use a lower tolerance for the q1 mesh because it always \n",
"# under-integrates due to the mesh 'cutting corners'\n",
"for index,mesh in enumerate([meshq1,meshq2]):\n",
" rtols = [1e-2,1e-5]\n",
" mesh.reset()\n",
" circleMesh(mesh)\n",
" checkSurfaceIntegration(mesh, 2.*math.pi, rtol=rtols[index]) \n",
" checkVolumeIntegration( mesh, math.pi, rtol=rtols[index])\n",
"\n",
"meshfig = vis.Figure(figsize=(600,600))\n",
"meshfig.append( vis.objects.Mesh(meshq1) )\n",
"meshfig.show()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\tGlobal element size: 8x8x8\n",
"\tLocal offset of rank 0: 0x0x0\n",
"\tLocal range of rank 0: 8x8x8\n",
"\tGlobal element size: 8x8x8\n",
"\tLocal offset of rank 0: 0x0x0\n",
"\tLocal range of rank 0: 8x8x8\n"
]
},
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# lets do some 3d tests\n",
"minCoord = [-5.,-2.,-1.]\n",
"maxCoord = [ 5., 2., 1.]\n",
"meshq1 = uw.mesh.FeMesh_Cartesian(elementType='Q1', elementRes=(8,8,8), minCoord=minCoord, maxCoord=maxCoord)\n",
"meshq2 = uw.mesh.FeMesh_Cartesian(elementType='Q2', elementRes=(8,8,8), minCoord=minCoord, maxCoord=maxCoord)\n",
"\n",
"def correctSurface(minCoord,maxCoord):\n",
" val = 0.\n",
" val += ( maxCoord[0] - minCoord[0] )*( maxCoord[1] - minCoord[1] ) \n",
" val += ( maxCoord[0] - minCoord[0] )*( maxCoord[2] - minCoord[2] ) \n",
" val += ( maxCoord[1] - minCoord[1] )*( maxCoord[2] - minCoord[2] ) \n",
" return 2*val\n",
"\n",
"for mesh in [meshq1,meshq2]:\n",
" mesh.reset()\n",
" checkSurfaceIntegration(mesh, correctSurface(minCoord,maxCoord)) \n",
" checkVolumeIntegration( mesh, correctVol(minCoord,maxCoord))\n",
"meshfig = vis.Figure(figsize=(800,400))\n",
"meshfig.append( vis.objects.Mesh(meshq1) )\n",
"meshfig.show()"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# check\n",
"for mesh in [meshq1,meshq2]:\n",
" mesh.reset()\n",
" rotate(mesh)\n",
" checkSurfaceIntegration(mesh, correctSurface(minCoord,maxCoord)) \n",
" checkVolumeIntegration( mesh, correctVol(minCoord,maxCoord))\n",
"meshfig.show()"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# refine, check, rotate, check\n",
"for mesh in [meshq1,meshq2]:\n",
" mesh.reset()\n",
" refine(mesh, minCoord, maxCoord)\n",
" checkSurfaceIntegration(mesh, correctSurface(minCoord,maxCoord)) \n",
" checkVolumeIntegration( mesh, correctVol(minCoord,maxCoord))\n",
" rotate(mesh)\n",
" checkSurfaceIntegration(mesh, correctSurface(minCoord,maxCoord)) \n",
" checkVolumeIntegration( mesh, correctVol(minCoord,maxCoord))\n",
"meshfig.show()"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\tGlobal element size: 16x16x16\n",
"\tLocal offset of rank 0: 0x0x0\n",
"\tLocal range of rank 0: 16x16x16\n",
"\tGlobal element size: 16x16x16\n",
"\tLocal offset of rank 0: 0x0x0\n",
"\tLocal range of rank 0: 16x16x16\n"
]
},
{
"data": {
"text/html": [
"
"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# giordani mesh \n",
"minCoord = [-1.,-1.,-1.]\n",
"maxCoord = [ 1., 1., 1.]\n",
"meshq1 = uw.mesh.FeMesh_Cartesian(elementType='Q1', elementRes=(16,16,16), minCoord=minCoord, maxCoord=maxCoord)\n",
"meshq2 = uw.mesh.FeMesh_Cartesian(elementType='Q2', elementRes=(16,16,16), minCoord=minCoord, maxCoord=maxCoord)\n",
"def sphereMesh(mesh):\n",
" with mesh.deform_mesh():\n",
" for vert in mesh.data[:]:\n",
" radius = np.max(np.abs(vert))\n",
" angle1 = np.arctan2(vert[1],vert[0])\n",
" angle2 = np.arctan2(vert[2],np.linalg.norm(vert[0:2]))\n",
" distinplane = radius*math.cos(angle2)\n",
" vert[0] = distinplane*math.cos(angle1)\n",
" vert[1] = distinplane*math.sin(angle1)\n",
" vert[2] = radius*math.sin(angle2)\n",
"\n",
"# check.. note we need to use a lower tolerance for the q1 mesh because it always \n",
"# under-integrates due to the mesh 'cutting corners'\n",
"for index,mesh in enumerate([meshq1,meshq2]):\n",
" rtols = [1e-2,1e-5]\n",
" mesh.reset()\n",
" sphereMesh(mesh)\n",
" checkSurfaceIntegration(mesh, 4.*math.pi, rtol=rtols[index]) \n",
" checkVolumeIntegration( mesh, 4./3.*math.pi, rtol=rtols[index])\n",
"\n",
"meshfig = vis.Figure(figsize=(600,600))\n",
"meshfig.append( vis.objects.Mesh(meshq1) )\n",
"meshfig.show()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# integrate radial function\n",
"import underworld.function as fn\n",
"coord = fn.input()\n",
"r = fn.math.sqrt( fn.math.dot(coord,coord) ) # exact solution for integration across unit sphere is \\pi\n",
"\n",
"for index,mesh in enumerate([meshq1,meshq2]):\n",
" rtols = [1e-2,1e-5]\n",
" checkVolumeIntegration( mesh, math.pi, fn=r, rtol=rtols[index])\n",
"\n",
"# lets also integrate just one eighth of sphere surface\n",
"func = fn.branching.conditional( [ ( (coord[0]>0.) & (coord[1]>0.) & (coord[2] > 0.), 1.),\n",
" ( True , 0.) ] ) \n",
"for index,mesh in enumerate([meshq1,meshq2]):\n",
" rtols = [1e-2,1e-5]\n",
" checkSurfaceIntegration( mesh, 1./2.*math.pi, fn=func, rtol=rtols[index])"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# finally let's add an empty surface integral tests.\n",
"# this is in particular important for parallel surface integrations\n",
"# where often some procs have empty sets. it should return\n",
"# a valid result, zero. \n",
"zero = uw.utils.Integral( fn=1., mesh=meshq1, integrationType='Surface', surfaceIndexSet=meshq1.specialSets[\"Empty\"])\n",
"result = zero.evaluate()\n",
"assert np.allclose( result, 0., rtol=1e-5 ), \"Error occurred in surface integration. Integration = {}, Expected = {}\".format(result[0], 0.)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.10.4"
}
},
"nbformat": 4,
"nbformat_minor": 4
}