3
$\begingroup$

enter image description hereenter image description hereenter image description hereNow I am try to convert matlab file to mathematica,

This is the original Matlab command, surf, for making surface plots:

surf (xei , yei , zei , EE,'FaceColor ','interp ','FaceAlpha ' ,1.0 , 'LineWidth ' ,... 301
0.01 , 'EdgeAlpha ' ,0.2 , 'EdgeColor ','w');

If I want to plot the same matrix in Mathematica, which correct Plot Function shall I use?

For a better understanding, I put some convert codes here.

n = 70;
m = n/2;
CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0, 
    0}, {63.9, 63.9, 165.7, 0, 0, 0}, {0, 0, 0, 2*79.6, 0, 0}, {0, 0, 
    0, 0, 2*79.6, 0}, {0, 0, 0, 0, 0, 2*79.6}};
phi = N[Range[0, Pi, Pi/(m - 1)]];
theta = N[Range[0, 2*Pi, 2*Pi/(n - 1)]];
EE = ConstantArray[0, {m, n}];

xei = ConstantArray[0, {m, n}];
yei = ConstantArray[0, {m, n}];
zei = ConstantArray[0, {m, n}];

SS = Inverse[CC];



For[i = 1, i < n + 1, i++, For[j = 1, j < m + 1, j++,
  xhen = {{Sin[phi[[j]]]*Cos[theta[[i]]], 
     Sin[phi[[j]]]*Sin[theta[[i]]], Cos[phi[[j]]]}};
  (*xheng=Transpose[xhen];*)
  xiaoa = KroneckerProduct[Transpose[xhen], xhen];
  dbV = {{xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]], 
     Sqrt[2]*xiaoa[[2, 3]], Sqrt[2]*xiaoa[[1, 3]], 
     Sqrt[2]*xiaoa[[1, 2]]}};
  (*dbV=Transpose[dbV];*)
  FFXX = SS . Transpose[dbV];
  (*FXXX=Transpose[dbV];*)
  EE[[j, i]] = 1/(dbV . FFXX);]]


For[i = 1, i < n + 1, i++, For[j = 1, j < m + 1, j++,
  xei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Cos[theta[[i]]];
  yei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Sin[theta[[i]]];
  zei[[j, i]] = EE[[j, i]]*Cos[phi[[j]]];]]
ListPlot3D[Transpose[{xei, yei, zei}]]
$\endgroup$
  • 1
    $\begingroup$ Try ListPlot3D[Transpose[{xei, yei, zei}]]. $\endgroup$ – Henrik Schumacher Jul 8 '18 at 13:25
  • $\begingroup$ @HenrikSchumacher, xei,yei,zei are all not vector/list, if I use your method, then ListPlot3D::arrayerr: {<<1>>} must be a valid array or a list of valid arrays. $\endgroup$ – ABCDEMMM Jul 8 '18 at 13:34
  • $\begingroup$ So, what is xei, yei, zei`? Example data is always welcome. $\endgroup$ – Henrik Schumacher Jul 8 '18 at 13:37
  • $\begingroup$ Hey @HenrikSchumacher, I put complete code here, thanks a lot for your support! $\endgroup$ – ABCDEMMM Jul 8 '18 at 13:47
  • 1
    $\begingroup$ Try also ListSurfacePlot3D (along with Flatten as in Anton Antonov's answer). $\endgroup$ – Henrik Schumacher Jul 8 '18 at 16:20
6
$\begingroup$

Edit 4

We can project it on xy plane by setting z=0.

n = 250;
m = n/2;
CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0, 
    0}, {63.9, 63.9, 165.7, 0, 0, 0}, {0, 0, 0, 2*79.6, 0, 0}, {0, 0, 
    0, 0, 2*79.6, 0}, {0, 0, 0, 0, 0, 2*79.6}};
phi = N[Range[0, Pi, Pi/(m - 1)]];
theta = N[Range[0, 2*Pi, 2*Pi/(n - 1)]];
EE = xei = yei = zei = ConstantArray[0, {m, n}];
SS = Inverse[CC];

Table[xhen = {Sin[phi[[j]]]*Cos[theta[[i]]], 
    Sin[phi[[j]]]*Sin[theta[[i]]], Cos[phi[[j]]]};
  xiaoa = KroneckerProduct[xhen, xhen];
  dbV = {xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]], 
    Sqrt[2]*xiaoa[[2, 3]], Sqrt[2]*xiaoa[[1, 3]], 
    Sqrt[2]*xiaoa[[1, 2]]};
  FFXX = SS.dbV;
  EE[[j, i]] = 1/(dbV.FFXX), {i, n}, {j, m}];

Table[xei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Cos[theta[[i]]];
  yei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Sin[theta[[i]]];
  zei[[j, i]] = EE[[j, i]]*Cos[phi[[j]]];, {i, n}, {j, m}];

data = Transpose[Flatten /@ {xei, yei, zei, EE}];
data[[All, 3]] = 0;
color = Blend["ThermometerColors", Rescale[#, MinMax@EE]] &;



Legended[Graphics3D[{color[Last[#]], AbsolutePointSize[20], 
     Point[Most[#]]} & /@ data, ImageSize -> Large, 
  ViewPoint -> Above, Boxed -> False], 
 BarLegend[{color, MinMax@EE}, LegendMarkerSize -> {22, 400}]]

enter image description here

Or

data = Transpose[Flatten /@ {xei, yei, zei, EE}];
data = data[[All, {1, 2, 4}]];

Legended[Graphics[{color[Last[#]], AbsolutePointSize[20], 
     Point[Most[#]]} & /@ data, ImageSize -> Large], 
 BarLegend[{color, MinMax@EE}, LegendMarkerSize -> {22, 400}]]

enter image description here

Or you achieve this by using @xzczd

CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0, 
    0}, {63.9, 63.9, 165.7, 0, 0, 0}, {0, 0, 0, 2 79.6, 0, 0}, {0, 0, 
    0, 0, 2 79.6, 0}, {0, 0, 0, 0, 0, 2 79.6}};
SS = Inverse[CC];
xhen = {Sin[phi] Cos[theta], Sin[phi] Sin[theta], Cos[phi]};
xiaoa = Outer[Times, xhen, xhen];
dbV = {xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]], 
   Sqrt[2] xiaoa[[2, 3]], Sqrt[2] xiaoa[[1, 3]], 
   Sqrt[2] xiaoa[[1, 2]]};
EE = 1/dbV.SS.dbV;

SphericalPlot3D[EE, {phi, 0, Pi}, {theta, 0, 2 Pi}, 
 ColorFunction -> 
  Function[{x, y, z, phi, theta, r}, 
   ColorData["ThermometerColors"][r]], ViewPoint -> Above, 
 Boxed -> False, Axes -> False, PlotPoints -> 100, Mesh -> None]

enter image description here

Edit 3

You can increase number of point by increasing n.

enter image description here

 {col, row} = ImageDimensions@parula;
    ParulaMMA = 
      Module[{colorlist}, 
       colorlist = 
        Catenate@
         ImageData@ImageTake[parula, {Round[row/2], Round[row/2]}, All];
       Evaluate[Blend[RGBColor @@@ colorlist, Rescale[#, MinMax@EE]] &]];



n = 250;
m = n/2;
CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0, 
    0}, {63.9, 63.9, 165.7, 0, 0, 0}, {0, 0, 0, 2*79.6, 0, 0}, {0, 0, 
    0, 0, 2*79.6, 0}, {0, 0, 0, 0, 0, 2*79.6}};
phi = N[Range[0, Pi, Pi/(m - 1)]];
theta = N[Range[0, 2*Pi, 2*Pi/(n - 1)]];
EE = xei = yei = zei = ConstantArray[0, {m, n}];
SS = Inverse[CC];

Table[
  xhen = {Sin[phi[[j]]]*Cos[theta[[i]]], 
    Sin[phi[[j]]]*Sin[theta[[i]]], Cos[phi[[j]]]};
  xiaoa = KroneckerProduct[xhen, xhen];
  dbV = {xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]], 
    Sqrt[2]*xiaoa[[2, 3]], Sqrt[2]*xiaoa[[1, 3]], 
    Sqrt[2]*xiaoa[[1, 2]]};
  FFXX = SS.dbV;
  EE[[j, i]] = 1/(dbV.FFXX), {i, n}, {j, m}];

Table[xei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Cos[theta[[i]]];
  yei[[j, i]] = EE[[j, i]]*Sin[phi[[j]]]*Sin[theta[[i]]];
  zei[[j, i]] = EE[[j, i]]*Cos[phi[[j]]];, {i, n}, {j, m}];

data = Transpose[Flatten /@ {xei, yei, zei, EE}];


Legended[Graphics3D[{ParulaMMA[Last[#]], AbsolutePointSize[20], 
     Point[Most[#]]} & /@ data, Axes -> True, ImageSize -> Large, 
  PlotRange -> {{-200, 200}, {-200, 200}, {-200, 200}}], 
 BarLegend[{ParulaMMA, MinMax@EE}, LegendMarkerSize -> {22, 400}]]

enter image description here

Or

Mathematica's built in gradient color

color = Blend["ThermometerColors", Rescale[#, MinMax@EE]] &;

Legended[Graphics3D[{color[Last[#]], AbsolutePointSize[20], 
     Point[Most[#]]} & /@ data, Axes -> True, ImageSize -> Large, 
  PlotRange -> {{-200, 200}, {-200, 200}, {-200, 200}}], 
 BarLegend[{color, MinMax@EE}, LegendMarkerSize -> {22, 400}]]

enter image description here

Edit 2 Here is another way to visualize it. see this for parula color.

data = Transpose[Flatten /@ {xei, yei, zei, EE}];

enter image description here

{col, row} = ImageDimensions@parula;
ParulaMMA = 
  Module[{colorlist}, 
   colorlist = 
    Catenate@
     ImageData@ImageTake[parula, {Round[row/2], Round[row/2]}, All];
   Evaluate[Blend[RGBColor @@@ colorlist, Rescale[#, MinMax@EE]] &]];

Legended[Graphics3D[{ParulaMMA[Last[#]], Sphere[Most[#], 20]} & /@ 
   data, Axes -> True, ImageSize -> Large, 
  PlotRange -> {{-200, 200}, {-200, 200}, {-200, 200}}], 
 BarLegend[{ParulaMMA, MinMax@EE}, LegendMarkerSize -> {22, 400}]]

enter image description here

Original Answer

This is the closest I can get.

data = Transpose[Flatten /@ {xei, yei, zei, EE}];

ListPointPlot3D[List /@ Most /@ data, 
 PlotStyle -> ({AbsolutePointSize[22], 
      Blend["BlueGreenYellow", Rescale[#, MinMax@EE]]} & /@ 
    Last /@ data), BoxRatios -> {1, 1, 1}]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks a lot for figuring this out! It looks much better, but surface is not smoothly in Mathematica. $\endgroup$ – ABCDEMMM Jul 9 '18 at 8:11
  • $\begingroup$ If I run this plot, I get bugs, namely: ImageDimensions::imginv: Expecting an image or graphics instead of parula. Set::shape: Lists {col,row} and ImageDimensions[parula] are not the same shape. ImageTake::imginv: Expecting an image or graphics instead of parula. $\endgroup$ – ABCDEMMM Jul 9 '18 at 15:09
  • $\begingroup$ You need to define parula. got to here mathworks.com/help/matlab/ref/colormap.html and copy parula and paste it in Mathematica notebook $\endgroup$ – Okkes Dulgerci Jul 9 '18 at 15:22
  • $\begingroup$ Great, thanks a lot, this problem is solved! $\endgroup$ – ABCDEMMM Jul 9 '18 at 15:43
  • $\begingroup$ @Okkes Dulgerci, is there also a simple command that we can get a 2D results(projection:3D-2D) $\endgroup$ – ABCDEMMM Aug 10 '18 at 11:05
5
$\begingroup$

I should have read your question and the surf help page more carefully -- your are using four matrices {xei, yei, zei, EE}.

Here is one way to use EE for the surface color:

colorRules = 
  Most[#] -> Last[#] & /@ Transpose[Flatten /@ {xei, yei, zei, EE}];

ListSurfacePlot3D[Transpose[Flatten /@ {xei, yei, zei}], 
 ColorFunction -> Function[{x, y, z}, {x, y, z} /. colorRules]]

enter image description here

I assume this can be improved.

Second try

I understand that the following plots are "volume" density and contour plots, not a surface density plot which you want. But I think they are close...

Using ListDensityPlot3D:

ListDensityPlot3D[Transpose[Flatten /@ {xei, yei, zei, EE}]]

enter image description here

Another plot that is close is made with ListContourPlot3D:

ListContourPlot3D[Transpose[Flatten /@ {xei, yei, zei, EE}], 
 Contours -> 30, Mesh -> None, PlotLegends -> "Expressions"]

enter image description here

First take

(More of a comment to find out what OP wants...)

Evaluating the code you posted we have these dimensions for {xei, yei, zei}:

Dimensions /@ {xei, yei, zei}    
(* {{35, 70, 1, 1}, {35, 70, 1, 1}, {35, 70, 1, 1}} *)

I strongly suspect that what you want is this plot:

ListPlot3D[Transpose[Flatten /@ {xei, yei, zei}], PlotRange -> All]

enter image description here

$\endgroup$
  • $\begingroup$ But It still contains the bugs as the cube is not closed form. looking for your reply! Thanks $\endgroup$ – ABCDEMMM Jul 8 '18 at 20:37
  • $\begingroup$ the mesh gitter also not correct. $\endgroup$ – ABCDEMMM Jul 8 '18 at 20:52
  • $\begingroup$ why 'ListSurfacePlot3D' is not correct in Mathematica? $\endgroup$ – ABCDEMMM Jul 9 '18 at 10:29
  • $\begingroup$ @FEAPMAN "why 'ListSurfacePlot3D' is not correct in Mathematica?" -- I do not know is it correct or not. Meaning, ListSurfacePlot3D might have a bug, or the problem might be in your data and a particular "feature" of ListSurfacePlot3D... It is a good idea to investigate, but I am not doing that anytime soon. $\endgroup$ – Anton Antonov Jul 9 '18 at 14:06
4
$\begingroup$

OK, though OP doesn't reply to my comment, I'd like to post a possible answer. If OP just wants to obtain that plot and "plot the same matrix" is not necessary, then the solution is quite simple. We don't need to discretize the analytic solution and transform the coordinate etc., because we have SphericalPlot3D:

CC = {{165.7, 63.9, 63.9, 0, 0, 0}, {63.9, 165.7, 63.9, 0, 0, 0}, {63.9, 63.9, 165.7, 0, 
    0, 0}, {0, 0, 0, 2 79.6, 0, 0}, {0, 0, 0, 0, 2 79.6, 0}, {0, 0, 0, 0, 0, 2 79.6}};
SS = Inverse[CC];
xhen = {Sin[phi] Cos[theta], Sin[phi] Sin[theta], Cos[phi]};
xiaoa = Outer[Times, xhen, xhen];
dbV = {xiaoa[[1, 1]], xiaoa[[2, 2]], xiaoa[[3, 3]], Sqrt[2] xiaoa[[2, 3]], 
   Sqrt[2] xiaoa[[1, 3]], Sqrt[2] xiaoa[[1, 2]]};
EE = 1/dbV.SS.dbV;

SphericalPlot3D[EE, {phi, 0, Pi}, {theta, 0, 2 Pi}, 
 ColorFunction -> Function[{x, y, z, phi, theta, r}, ColorData["AvocadoColors"][r]]]

Mathematica graphics

I believe the definition of EE can be simplified further, but without background information it's somewhat beyond my reach.

$\endgroup$
  • $\begingroup$ this plot results can show us the Anisotropic Elasticity (similar to the work: "Anisotropic Elasticity" demonstrations.wolfram.com/AnisotropicElasticity Wolfram Demonstrations Project Published: March 7 2011) $\endgroup$ – ABCDEMMM Jul 10 '18 at 12:10
  • $\begingroup$ The original code is wrotten in Matlab which needs to discretize the analytic solution and transform the coordinate etc. Now Mathematica show its strongly numerical ability in this case, cheers! $\endgroup$ – ABCDEMMM Jul 10 '18 at 12:13
  • $\begingroup$ @FEAPMAN Is EE defined by some kind of rotation? Can you elaborate a little on this part? $\endgroup$ – xzczd Jul 10 '18 at 15:05
  • $\begingroup$ EE is elasticity tensor. $\endgroup$ – ABCDEMMM Aug 10 '18 at 11:03
  • $\begingroup$ @ABCDEMMM I know, but what I'm interested in is, why is it defined in this way?xiaoa and dbV looks like something related to rotation or coordinate transform, and I believe they can be simplified, but without more background information it's beyond my reach. $\endgroup$ – xzczd Aug 11 '18 at 14:21

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.