3
$\begingroup$

I have made a code for the solution to a particular problem involving the Laplacian, and after plotting in 3D the resultant potential, I have tried to look for a way to extract points from that graph in order to later get the electric field produce by that potential, so I wonder if there is a way to do this, or another way to plot the electric field by getting the following solution to the Laplacian?

This is the solution

a = 4;
f[x_, y_] := Sin[ArcTan[Abs[y/x]]]
g[x_, y_] := -Sin[ArcTan[Abs[y/x]]]
leqn = Laplacian[u[x, y], {x, y}] == 0;
\[CapitalOmega] = 
  ImplicitRegion[x^2 + y^2 >= 1, {{x, -a, a}, {y, -a, a}}];
A[x_, y_] := NDSolveValue[{\!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(u[x, y]\)\) == 0, 
   DirichletCondition[u[x, y] == f[x, y], y == Sqrt[1 - x^2]], 
   DirichletCondition[u[x, y] == g[x, y], y == -Sqrt[1 - x^2]], 
   DirichletCondition[u[x, y] == 0, 
    x == a || x == -a || y == a || y == -a]}, 
  u, {x, y} \[Element] \[CapitalOmega]]
ContourPlot[A[x, y], {x, y} \[Element] \[CapitalOmega], 
 PlotRange -> All, FrameLabel -> {x, y}, PlotPoints -> 20, 
 ColorFunction -> "Rainbow", PlotLegends -> Automatic]
Plot3D[A[x, y], {x, y} \[Element] \[CapitalOmega], PlotRange -> All, 
 ColorFunction -> "Rainbow", ViewPoint -> {0, 0, \[Infinity]}, 
 PlotLegends -> Automatic]

The ContourPlot

The Plot3D

So I tried using this idea from http://community.wolfram.com/groups/-/m/t/452396, by putting Array[{##, A[##]} &, {10, 10}, {x, y} \[Element] \[CapitalOmega]] Graphics3D[Point /@ %]

But it doesn't seem to work, or I am using it wrong. Does anyone has an idea of how to solve this problem?

$\endgroup$
3
  • $\begingroup$ Does Array[{##, A[##]} &, {10, 10}, {{-a, a}, {-a, a}}]; Graphics3D[Point /@ %] give what you need? $\endgroup$ – kglr Oct 25 '17 at 6:31
  • $\begingroup$ It doesn't seem to work, it gives me two errors: InterpolatingFunction::femdmval: Input value {-0.444444,-0.444444} lies outside the range of data in the interpolating function. And, General::stop: Further output of InterpolatingFunction::femdmval will be suppressed during this calculation. $\endgroup$ – Mounice Oct 25 '17 at 6:34
  • $\begingroup$ Are you sure about that Sin[ArcTan[Abs[y/x]]? It seems to me you should be using two-argument arctangent. $\endgroup$ – J. M.'s torpor Oct 25 '17 at 7:12
3
$\begingroup$
aa = NDSolveValue[{leqn, 
     DirichletCondition[u[x, y] == f[x, y], y == Sqrt[1 - x^2]], 
     DirichletCondition[u[x, y] == g[x, y], y == -Sqrt[1 - x^2]], 
     DirichletCondition[u[x, y] == 0, x == a || x == -a || y == a || y == -a]}, 
   u, Element[{x, y}, Ω]];

points = Array[If[RegionMember[Ω, {##}], {##, aa[##]} , ##&[]]&, 
 {10, 10}, {{-a, a}, {-a, a}}];
Graphics3D[{Red, PointSize[Large], Point /@ points}, BoxRatios -> 1]

enter image description here

Show[Plot3D[aa[x, y], Element[{x, y}, Ω],  PlotRange -> All, 
   ColorFunction -> "Rainbow",  PlotLegends -> Automatic], g3d]

enter image description here

$\endgroup$
4
  • $\begingroup$ Using RegionMember in this way can be very expensive. You could generate a rmf=RegionMember[\[CapitalOmega]] and use that in the If statement. But using a discretized version of Omega would be even better. $\endgroup$ – user21 Oct 28 '17 at 0:00
  • $\begingroup$ @user21, thank you. I will update with your suggestion. $\endgroup$ – kglr Oct 28 '17 at 13:56
  • $\begingroup$ @kglr is there a new command in MMA 12.2 for this task? $\endgroup$ – ABCDEMMM Jun 7 at 1:01
  • $\begingroup$ @ABCDEMMM, i don't know if there is a new command that does this. $\endgroup$ – kglr Jun 7 at 9:21
3
$\begingroup$

Here is another alternative:

A = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0, 
      DirichletCondition[u[x, y] == f[x, y], y == Sqrt[1 - x^2]], 
      DirichletCondition[u[x, y] == g[x, y], y == -Sqrt[1 - x^2]], 
      DirichletCondition[u[x, y] == 0, 
        x == a || x == -a || y == a || y == -a]}, 
    u, {x, y} \[Element] \[CapitalOmega]]

We extract the values and the mesh (and from that the coordinates) from the interpolating function:

vals = A["ValuesOnGrid"];
mesh = A["ElementMesh"];
coords = mesh["Coordinates"];
data = Join[coords, Partition[vals, 1], 2];
Graphics3D[{Red, Point[data]}]

enter image description here

$\endgroup$

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.