This question is partly based off Gonza's question and Mauricio's answer. I wanted to adapt the FEM method to simulate a point pressure being applied to distort an encapsulated material that is resting on a surface. I'll include an image of what I've achieved so far, as I'm probably not describing it very well:
I'm trying to apply different material properties to both regions of the generated mesh, following roughly the example listed in ElementMeshCreation here. The code uses a replacement rule so that different properties are called depending on the value of ElementMarker with:
(*If statement to replace material properties depending on ElementMarker*)
props = If[ElementMarker == 10, Propiedades10, Propiedades20];
(*NDSolve with boundary conditions included*)
{u1, v1, \[Sigma]x1, \[Sigma]y1, \[Tau]xy1} = NDSolveValue[{PS == bcN, hl, bccorrect} /. props, {u, v, \[Sigma]x, \[Sigma]y, \[Sigma]xy}, Element[{x, y}, Mesh4]];
but I am getting this error:
If[ElementMarker==10,Propiedades10,Propiedades20]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing
Is it the wrong place to be inserting a replacement rule? I feel that the if statement may not be referencing the correct RegionMarker value? The full code I've used/modified is:
(*Initiation*)
Needs["NDSolve`FEM`"];
(*Forces*)
q = 60;
(*Material properties*)
Propiedades10 = {Y -> 25000000, \[Nu] -> 47/100};
Propiedades20 = {Y -> 105000000000, \[Nu] -> 30/100};
(*Geometry*)
(*https://reference.wolfram.com/language/FEMDocumentation/tutorial/\
ElementMeshCreation.html*)
L = 3;h1 = 0.4;sxMin = 0.1 L; sxMax = 0.9 L; syMin = 2 h1/5; syMax = 3 h1/5;
(*Mesh*)
Reg4 = ToBoundaryMesh["Coordinates" -> {{0., 0.}, {L, 0.}, {L, h1}, {0, h1}, {sxMin,syMin}, {sxMax, syMin}, {sxMax, syMax}, {sxMin, syMax}},"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4,1}, {5, 6}, {6, 7}, {7, 8}, {8, 5}}]}];
Mesh4 = ToElementMesh[Reg4,"RegionMarker" -> {{{L/2, h1/2}, 10},{{sxMin/2, syMin/2}, 20}},MeshQualityGoal -> 0];
(*FEM Solution*)
(*2D Hooke's law*)
hl = {\[Sigma]x[x, y] ==Y/(1 - \[Nu]^2) (D[u[x, y], x] + \[Nu] D[v[x, y], y]), \[Sigma]y[x, y] == Y/(1 - \[Nu]^2) (D[v[x, y], y] + \[Nu] D[u[x, y], x]), \[Sigma]xy[x, y] == (Y*\[Nu])/(1 - \[Nu]^2) (D[u[x, y], x] + D[v[x, y], y])};
(*Equations*)
PS = {Inactive[Div][{{0, -((Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + Inactive[Div][{{-(Y/(1 - \[Nu]^2)), 0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}.Inactive[Grad][u[x, y], {x, y}], {x, y}], Inactive[Div][{{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((Y \[Nu])/(1 - \[Nu]^2)),0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] + Inactive[Div][{{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][v[x, y], {x, y}], {x, y}]};
(*BCs*)
(*Neumann - Boundary condition of force exerted*)
bcN = {0, NeumannValue[-q, {1.25 <= x <= 1.75, y == h1}]};
(*Dirichlet - Boundary condition of held positions*)
bccorrect = {DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x <= 3 && y <= 0.01],DirichletCondition[v[x, y] == 0, x == L && y <= 0.01]};
(*If statement to replace material properties depending on ElementMarker*)
props = If[ElementMarker == 10, Propiedades10, Propiedades20];
(*NDSolve with boundary conditions included*)
{u1, v1, \[Sigma]x1, \[Sigma]y1, \[Tau]xy1} = NDSolveValue[{PS == bcN, hl, bccorrect} /. props, {u, v, \[Sigma]x, \[Sigma]y, \[Sigma]xy}, Element[{x, y}, Mesh4]];
(*Deformation*)
DMesh1 = ElementMeshDeformation[Mesh4, {u1, v1}, "ScalingFactor" -> 6*10^4];
(*Show deformed mesh in red and green*)
Show[{DMesh1["Wireframe"["MeshElementStyle" ->{Directive[FaceForm[Green]],Directive[FaceForm[Red]]}]]}, ImageSize -> Large]
Thanks!
*Edited code for clarity
*Edit 2: Line breaks in copy-pasted code were breaking it. Fixed this issue.
*Edit 3: Question wasn't clear. Edited.
