2
$\begingroup$

I have started to work on the examples presented in the SOLID MECHANICS technical notes. The example I am focusing on is the small, cantilever beam, with a load applied to the end surface.

enter image description here

I have modified various parameters from the original example to more closely resemble slightly different loading, material properties, etc..

enter image description here

Now that I am getting a handle on this - I want to see the mesh that is used in the solution process and get some statistics of the mesh used - number of elements, number of nodes, type of elements, etc.. Neither the first cantilever beam or the bracket Solid Mechanics example show how to plot out the mesh used in the solution process. This is probably an easy solution, but I just do not see it.

I have tried some approaches to extract the mesh - using the Mathematica Notebook Assistant, but not have any success. I have copied the code below - including some approaches I tried but have been commented out as they did not work. I do get a mesh plotted, but this mesh is "not valid" ?

enter image description here

Needs["NDSolve`FEM`"] $HistoryLength = 0; \[CapitalOmega] = Cuboid[{0, 0, 0}, {10, 1, 1}]; vars = {{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}}; (*pars=<|"Material"->Entity["Element","Titanium"]|>;*) pars = <|"Material" -> Entity["Element", "Aluminium"]|>; (*Subscript[\[CapitalGamma],force]=SolidBoundaryLoadValue[x==10,vars,\ pars,<|"Force"->{0,0,Quantity[-1000,"Newtons"]}|>];*) (*Subscript[\[CapitalGamma],force]=SolidBoundaryLoadValue[x==10,vars,\ pars,<|"Force"->{0,Quantity[-1000,"Newtons"],0}|>];*) (*Subscript[\[CapitalGamma],force]=SolidBoundaryLoadValue[z==1,vars,\ pars,<|"Force"->{0,0,Quantity[-1000,"Newtons"]}|>];*) (*Subscript[\[CapitalGamma],force]=SolidBoundaryLoadValue[z==1,vars,\ pars,<|"Pressure"->{0,0,Quantity[-1000,"Pascals"]}|>];*) Subscript[\[CapitalGamma], force] = SolidBoundaryLoadValue[boundaryLoadPosition, vars, pars, <|"Pressure" -> {0, Quantity[1000, "Pascals"], 0}|>]; (* y direction - shearing force *) (*boundaryLoadPosition=x>=0.01&&z==1;*) (*boundaryLoadPosition=0<=x<=10&&z==1;*) boundaryLoadPosition = 5 <= x <= 10 && 0.5 <= y <= 1 && z == 1; (*Subscript[\[CapitalGamma],force]=SolidBoundaryLoadValue[\ boundaryLoadPosition,vars,pars,<|"Pressure"->{0,0,Quantity[-1000,\ "Pascals"]}|>]; *)(* z direction*) backconstraintPositionsak = x == 0; Subscript[\[CapitalGamma], wall] = SolidFixedCondition[backconstraintPositionsak, vars, pars]; (* outline=Show[edgeframe,bmesh["Wireframe"["MeshElementStyle"\[Rule] Directive[Opacity[0.2],FaceForm[LightBlue],EdgeForm[]]]]]*) Graphics3D[{EdgeForm[], FaceForm[LightBlue], Opacity[0.2], Cuboid[{0, 0, 0}, {10, 1, 1}]}] outlinesak = Show[Graphics3D[{EdgeForm[], FaceForm[LightBlue], Opacity[0.2], \[CapitalOmega]}]]; outlinesakload = Show[Graphics3D[{EdgeForm[], FaceForm[LightBlue], Opacity[0.2], \[CapitalOmega]}, PlotLabel -> "Loading Area in RED"]]; outlinesakconstraint = Show[Graphics3D[{EdgeForm[], FaceForm[LightBlue], Opacity[0.2], \[CapitalOmega]}, PlotLabel -> "Constraint Area in BLUE"]]; (*boundaryLoadPosition=x>=0.01&&z==1;*) Show[outlinesakload, Graphics3D[{{RGBColor[1, 0, 0], Arrow[{{0, 0, 0}, {0.04`, 0, 0}}]}, {RGBColor[0, 1, 0], Arrow[{{0, 0, 0}, {0, 0.04`, 0}}]}, {RGBColor[0, 0, 1], Arrow[{{0, 0, 0}, {0, 0, 0.04`}}]}, {Text["X", {0.04`, 0, 0}], Text["Y", {0, 0.04`, 0}], Text["Z", {0, 0, 0.04`}]}}], Region[Style[ ImplicitRegion[ boundaryLoadPosition, {{x, -0.01, 10}, {y, 0, 1}, {z, 0, 1}}], Red]]] (*back contraint position*) Show[outlinesakconstraint, Graphics3D[{{RGBColor[1, 0, 0], Arrow[{{0, 0, 0}, {0.04`, 0, 0}}]}, {RGBColor[0, 1, 0], Arrow[{{0, 0, 0}, {0, 0.04`, 0}}]}, {RGBColor[0, 0, 1], Arrow[{{0, 0, 0}, {0, 0, 0.04`}}]}, {Text["X", {0.04`, 0, 0}], Text["Y", {0, 0.04`, 0}], Text["Z", {0, 0, 0.04`}]}}], Region[Style[ ImplicitRegion[ backconstraintPositionsak, {{x, 0, 10}, {y, 0, 1}, {z, 0, 1}}], Blue]]] MatrixForm[op = SolidMechanicsPDEComponent[vars, pars]] beamDisplacement = NDSolveValue[{op == Subscript[\[CapitalGamma], force], Subscript[\[CapitalGamma], wall]}, {u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z} \[Element] \[CapitalOmega]]; VectorDisplacementPlot3D[beamDisplacement, {x, y, z} \[Element] \[CapitalOmega], AxesLabel -> {"x", "y", "z"}, PlotLabel -> "Displacement Plot (units: XXX)"] VectorDisplacementPlot3D[ beamDisplacement, {x, y, z} \[Element] \[CapitalOmega], AxesLabel -> {"x", "y", "z"}, PlotLabel -> "Displacement Plot (units: XXX)", ViewVector -> {{-10, 0, 0}, {0, 0, 0}} (* Viewing from the negative x-axis towards the origin *) ] VectorDisplacementPlot3D[ beamDisplacement, {x, y, z} \[Element] \[CapitalOmega], AxesLabel -> {"x", "y", "z"}, PlotLabel -> "Displacement Plot (units: XXX)", ViewVector -> {{0, -10, 0}, {0, 0, 0}} (* Viewing from the negative y-axis towards the origin *) ] VectorDisplacementPlot3D[ beamDisplacement, {x, y, z} \[Element] \[CapitalOmega], AxesLabel -> {"x", "y", "z"}, PlotLabel -> "Displacement Plot Perspective (units: XXX)", ViewVector -> {{0, 0, -10}, {0, 0, 0}}, ViewVertical -> {0, 1, 0} (* Sets the y-axis as vertical *) ] VectorDisplacementPlot3D[ beamDisplacement, {x, y, z} \[Element] \[CapitalOmega], AxesLabel -> {"x", "y", "z"}, PlotLabel -> "Displacement Plot Orthographic (units: XXX)", ViewVector -> {{0, 0, -10}, {0, 0, 0}}, ViewVertical -> {0, 1, 0}, ViewProjection -> "Orthographic" (* Ensure orthographic projection is used *) ] (* Solve the FEM problem and extract the mesh *) (*beamDisplacementmesh=NDSolveValue[ { SolidMechanicsPDEComponent[vars,pars]==Subscript[\[CapitalGamma],\ force], Subscript[\[CapitalGamma],wall] }, {u[x,y,z],v[x,y,z],w[x,y,z]}, {x,y,z} \[Element] \[CapitalOmega], Method->{"PDEDiscretization"->{"FiniteElement", "MeshOptions"->{"MaxCellMeasure"->0.1}}} (* \ Control mesh refinement *) ];*) (*(* Extract the mesh from the beamDisplacement solution *) mesh=beamDisplacement["ElementMesh"]; (* Visualize the mesh *) Show[MeshRegion[mesh],MeshStyle->Red]*) (*(*Solve the FEM problem with mesh refinement*)beamDisplacementmesh=\ NDSolveValue[{SolidMechanicsPDEComponent[vars,pars]==Subscript[\ \[CapitalGamma],force],Subscript[\[CapitalGamma],wall]},{u[x,y,z],v[x,\ y,z],w[x,y,z]},{x,y,z}\[Element]\[CapitalOmega],Method->\ {"PDEDiscretization"->{"FiniteElement","MeshOptions"->\ {"MaxCellMeasure"->0.1}}}]; (*Extract the mesh from the solution*) mesh=beamDisplacementmesh["ElementMesh"]; (*Verify the mesh*) If[Head[mesh]===ElementMesh,(*Visualize the mesh*)Show[MeshRegion[\ mesh,MeshCellStyle->{{3,All}->FaceForm[Opacity[0.2,Red]],{2,All}->\ EdgeForm[Black]}],Axes->True,AxesLabel->{"x","y","z"},PlotLabel->\ "Finite Element Mesh"],Print["Error: Mesh extraction failed. Check \ the solution."]]*) meshsak = ToElementMesh[\[CapitalOmega]]["Wireframe"]; Show[meshsak] (* Number of elements *) numElements = MeshElementCount[meshsak] (* Type of elements *) elementTypes = MeshElementType[meshsak] (* Number of nodes *) numNodes = Length[MeshCoordinates[meshsak]] (* Check if meshsak is a valid ElementMesh *) If[Head[meshsak] === ElementMesh, Module[{numElements, elementTypes, numNodes}, (* Number of elements *) numElements = MeshElementCount[meshsak]; (* Type of elements *) elementTypes = MeshElementType[meshsak]; (* Number of nodes *) numNodes = Length[MeshCoordinates[meshsak]]; (* Display the statistics *) {numElements, elementTypes, numNodes} ], Print["Error: meshsak is not a valid ElementMesh."] ]
Improve this question
$\endgroup$

1 Answer 1

Reset to default
2
$\begingroup$

So, even the LLM created by WRI generates fake code ¯\_(ツ)_/¯. Just check document of ElementMesh and ToElementMesh by pressing F1, you'll find all you need… OK, this may not be immediately obvious, so let me post an answer. The Mathematica Notebook Assistant generated code involves the following mistakes:

  1. "Wireframe" is a "method" of ElementMesh[…] for visualization that generates a Graphics[…] or Graphics3D[…], so meshsak = ToElementMesh[Ω]["Wireframe"] should be meshsak = ToElementMesh[Ω]. See also the code in your previous questions for the correct usage of "Wireframe".

  2. There's no MeshElementCount function.

  3. MeshElementType is an option of ToElementMesh, not a function.

  4. MeshCoordinates doesn't have effect on ElementMesh[…].

So the correct code should be:

meshsak = ToElementMesh[Ω]; meshsak["Wireframe"] (* Graphics3D[…] *) meshsak@"NumberOfMeshElements" (* 9300 *) Head /@ meshsak["MeshElements"] (* {HexahedronElement} *) meshsak@"NumberOfCoordinates" (* 43307 *)

BTW, by checking

meshsak@"Methods"

you'll see all the available methods of ElementMesh[…].

As to the code in the comment, you just need to know a bit more about the list manipulation in Mathematica:

beamDisplacementmesh[[1, 0]]@"ElementMesh" (* ElementMesh[…] *) (* Alternatively: *) Head[First@beamDisplacementmesh]@"ElementMesh" (* ElementMesh[…] *)

You may read e.g. this and this post for more info. You should find more by searching in this site.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.