4
$\begingroup$

Modelling a Poisson Equation on piecewise constant material, ElementMarkers on the PointElements are used for the boundary conditions (using InitializeBoundaryConditions). This works fine.

In the same way I want to use markers on the MeshElements to determine the DiffusionCoefficient matrix on the specific element.

mesh with marker=1 within Rectangle[{-0.5,-0.5},{0.5,0.5}]

How can one retrieve the ElementMarkers in InitializePDECoefficients to implement a piecewise constant diffusion coefficient?

Based on the implicit region definition's inequalities, this code works:

pdeCoefficients = InitializePDECoefficients[
          variableData, solutionData, 
          "DiffusionCoefficients" -> {{
              If[-0.5 <= x <= 0.5 && -0.5 <= y <= 0.5, 
                   {{-11, 0}, {0, -11}}, {{-1, 0}, {0, -1}}] }}, 
          "DampingCoefficients" -> {{1}}]

The following code does run but the result is not the same. The diffusion matrix seems to be equal to the unity matrix everywhere.

sigma = Which[
   ElementMarker == 0, {{-1, 0}, {0, -1}}, 
   ElementMarker == 1, {{-11, 0}, {0, -11}}]; 
pdeCoefficients = 
   InitializePDECoefficients[
      variableData, solutionData, 
            "DiffusionCoefficients" -> {{sigma}}, 
            "DampingCoefficients" -> {{1}}]

What syntax options are allowed in DiffusionCoefficients? Does anybody have a description on the FEM package that is more in-depth than Wolfram Reference?

|improve this question
$\endgroup$
  • $\begingroup$ There is an example in the Markers section in the ElementMesh generation tutorial that shows that. If that does not help you'd need to clarify what exactly the issue is you are having by giving a complete code. $\endgroup$ – user21 Mar 24 '16 at 22:35
  • $\begingroup$ Where you able to figure this out? $\endgroup$ – user21 Mar 29 '16 at 8:05
  • $\begingroup$ I have already been through the Marker section of the ElementMesh tutorial. This does not really help. The cause of the trouble could be identified in the meantime: 1: MeshElement, BoundaryElement and PointElement markers do not operate independently (as stipulated in the 3D example of reference.wolfram.com/language/FEMDocumentation/tutorial/…). So all marker values must be disjoint where 0 must be assigned to "nothing". 2: The marker value "1" is automatically assigned to BoundaryElements denoting outer boundaries. This is not documented anywhere... $\endgroup$ – Harald F. Merkel Mar 29 '16 at 9:38
  • $\begingroup$ Could you tell me which paragraph does stipulate this for you. Then I could try to improve that. Also I am not sure I can follow you when you say that markers must be disjoint. There is no automatic marker assignment on boundaries except during the conversion from the boundary mesh to the full mesh. Also, I am not sure I understand how that relates to the question you posted here. If you could help understand that a bit better I could improve the documentation for future use if necessary. $\endgroup$ – user21 Mar 29 '16 at 21:06
3
$\begingroup$

Here is how you can do that:

Needs["NDSolve`FEM`"]
bmesh = ToBoundaryMesh[
   "Coordinates" -> {{0., 0.}, {1., 0.}, {1., 1.}, {0., 1.}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 1}, {3, 
        4}, {4, 1}}]}];
mesh = ToElementMesh[bmesh, 
   "RegionMarker" -> {{{1/2, 1/4}, 2}, {{1/2, 3/4}, 1}}, 
   MaxCellMeasure -> 0.05];
Show[
 mesh["Wireframe"],
 mesh["Wireframe"["MeshElementMarkerStyle" -> Blue]]
 ]

enter image description here

And then:

nr = ToNumericalRegion[mesh];
vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData["Space" -> nr];
sigma = Which[ElementMarker == 1, {{-1, 0}, {0, -1}}, 
   ElementMarker == 2, {{-11, 0}, {0, -11}}];
ipde = InitializePDECoefficients[vd, sd, 
   "DiffusionCoefficients" -> {{sigma}}, 
   "LoadCoefficients" -> {{1}}];
ibcs = InitializeBoundaryConditions[vd, 
   sd, {{DirichletCondition[u[x, y] == 0., y == 0]}}];
md = InitializePDEMethodData[vd, sd];
dpde = DiscretizePDE[ipde, md, sd];
dbcs = DiscretizeBoundaryConditions[ibcs, md, sd];
{l, s, d, m} = dpde["SystemMatrices"];
DeployBoundaryConditions[{l, s}, dbcs]
res = LinearSolve[s, l];
ifun = ElementMeshInterpolation[{mesh}, res];
Plot3D[ifun[x, y], {x, y} \[Element] mesh]

enter image description here

Looks good to me.

Update

In summary DirichletCondition boundary look at markers in PointElements, NeumanValue and PeriodicBoundaryCondition look at markers in Boundary Elements and the PDE coeffcients look at markers in mesh elements.

|improve this answer
$\endgroup$
  • $\begingroup$ So in summary Dirichlet boundary look at PointElement markers, Neumanvalues look at Boundary Elements markers and the PDE coeffcients look at RegionMarkers? That wasn't obvious at all to me. $\endgroup$ – Fortsaint Dec 16 '16 at 6:09
  • $\begingroup$ @Fortsaint, yes that's correct. I added a section to the documentation of the Element Mesh Generation that explains that better and will become available in the next release (V11.1). I hope that clarifies things. If you find other FEM stuff that's not clear please don't hesitate to ping me about it and I'll see what I can do. $\endgroup$ – user21 Dec 16 '16 at 10:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.