3D_Integrals

{  3D_INTEGRALS.PDE

This problem demonstrates the specification of various integrals in 3D.

( This is a modification of problem  3D_BRICKS.PDE)

}

title '3D Integrals'

coordinates

   cartesian3

variables

   Tp

definitions

   long = 1

   wide = 1

   K                           { thermal conductivity -- values supplied later }

   Q = 10*max(1-x^2-y^2-z^2,0)              { Thermal source }

   flag22=0        { build a test function for region 2, layer 2 }

   check22 = if flag22>0 then Tp else 0

   flag20=0        { build a test function for region 2, all layers }

   check20 = if flag20>0 then Tp else 0

initial values

   Tp = 0.

equations

   div[k*grad(Tp)] + Q = 0        { the heat equation }

extrusion

surface "bottom" z = -long

layer 'bottom'

surface "middle" z=0

layer 'top'

surface 'top' z= long        { divide Z into two layers }

boundaries

   surface 1 value(Tp)=0        { fix bottom surface temp }

   surface 3 value(Tp)=0        { fix top surface temp }

   Region 1                { define full domain boundary in base plane }

      layer 1 k=1                { bottom right brick }

      layer 2 k=0.1                { top right brick }

      start "outside" (-wide,-wide)

        value(Tp) = 0                { fix all side temps }

        line to (wide,-wide)        { walk outer boundary in base plane }

          to (wide,wide)

          to (-wide,wide)

          to close

   Region 2 "Left"                { overlay a second region in left half }

      flag20=1

      layer 1 k=0.2                { bottom left brick }

      layer 2 k=0.4  flag22=1        { top left brick }

      start(-wide,-wide)

        line to (0,-wide)                { walk left half boundary in base plane }

          to (0,wide)

          to (-wide,wide)

          to close

monitors

     contour(Tp) on surface z=0  as "XY Temp"

     contour(Tp) on surface x=0  as "YZ Temp"

     contour(Tp) on surface y=0  as "ZX Temp"

     elevation(Tp) from (-wide,0,0) to (wide,0,0)  as "X-Axis Temp"

     elevation(Tp) from (0,-wide,0) to (0,wide,0)  as "Y-Axis Temp"

     elevation(Tp) from (0,0,-long) to (0,0,long)  as "Z-Axis Temp"

plots

     contour(Tp) on z=0  as "XY Temp"

     contour(Tp) on x=0  as "YZ Temp"

     contour(Tp) on y=0  as "ZX Temp"

     contour(k*dz(Tp)) on z=-0.001  as "Low Middle Z-Flux"

     contour(k*dz(Tp)) on z=0.001  as "High Middle Z-Flux"

     summary

                report(integral(Tp,2,2))                { integrals in region 2, layer 2 }

       report(integral(Tp,"Left","Top"))        

       report(integral(check22))

       report '-----'

                report(integral(Tp,2,0))                { integrals in region 2, all layers (selected by 0) }

                report(integral(check20))

       report '-----'

                report(integral(Tp,"ALL","ALL"))        { integrals over total volume }

                report(integral(Tp))

       report '-----'

       report(sintegral(normal(-k*grad(Tp)),2))        { surface integral on 'middle' }

       report(sintegral(normal(-k*grad(Tp)),'Middle'))        { surface integral on 'middle' }

       report(sintegral(-k*dz(Tp),2))        { surface integral on 'middle' }

       report '-----'

       report(sintegral(normal(-k*grad(Tp)),1))        as "Bottom Flux"

       report(sintegral(normal(-k*grad(Tp)),3))        as "Top Flux"

       report(sintegral(normal(-k*grad(Tp)),"outside")) as "Side Flux"

       report(sintegral(normal(-k*grad(Tp)),1)+sintegral(normal(-k*grad(Tp)),3)

               +sintegral(normal(-k*grad(Tp)),"outside")) as "Bottom+Top+Side Flux"

       report(sintegral(normal(-k*grad(Tp))))        { surface integral on total outer surface }

       report(integral(Q) )        as "Source Integral"

       report(sintegral(normal(-k*grad(Tp)),"outside")

               +sintegral(normal(-k*grad(Tp)),1)+sintegral(normal(-k*grad(Tp)),3)

               -integral(Q) )        as "Energy Error"

       report(integral(div(-k*grad(Tp))) )        as "Divergence Integral"

       report '-----'

       { surface integral over outer surface of region 2, layer 2 }

       report(sintegral(normal(-k*grad(Tp)),"Left","Top"))        

       report(integral(Q,"Left","Top"))        

       report '-----'

end