5
$\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}]]
Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ Try ListPlot3D[Transpose[{xei, yei, zei}]]. $\endgroup$
    – Henrik Schumacher
    Jul 8, 2018 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, 2018 at 13:34
  • $\begingroup$ So, what is xei, yei, zei`? Example data is always welcome. $\endgroup$
    – Henrik Schumacher
    Jul 8, 2018 at 13:37
  • $\begingroup$ Hey @HenrikSchumacher, I put complete code here, thanks a lot for your support! $\endgroup$
    – ABCDEMMM
    Jul 8, 2018 at 13:47
  • 1
    $\begingroup$ Try also ListSurfacePlot3D (along with Flatten as in Anton Antonov's answer). $\endgroup$
    – Henrik Schumacher
    Jul 8, 2018 at 16:20

3 Answers 3

Reset to default
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

Improve this answer
$\endgroup$
14
  • $\begingroup$ Thanks a lot for figuring this out! It looks much better, but surface is not smoothly in Mathematica. $\endgroup$
    – ABCDEMMM
    Jul 9, 2018 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, 2018 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$
    – OkkesDulgerci
    Jul 9, 2018 at 15:22
  • $\begingroup$ Great, thanks a lot, this problem is solved! $\endgroup$
    – ABCDEMMM
    Jul 9, 2018 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, 2018 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

Improve this answer
$\endgroup$
4
  • $\begingroup$ But It still contains the bugs as the cube is not closed form. looking for your reply! Thanks $\endgroup$
    – ABCDEMMM
    Jul 8, 2018 at 20:37
  • $\begingroup$ the mesh gitter also not correct. $\endgroup$
    – ABCDEMMM
    Jul 8, 2018 at 20:52
  • $\begingroup$ why 'ListSurfacePlot3D' is not correct in Mathematica? $\endgroup$
    – ABCDEMMM
    Jul 9, 2018 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, 2018 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.

Improve this answer
$\endgroup$
5
  • $\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, 2018 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, 2018 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, 2018 at 15:05
  • $\begingroup$ EE is elasticity tensor. $\endgroup$
    – ABCDEMMM
    Aug 10, 2018 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, 2018 at 14:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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