10
$\begingroup$

I simulated some nonlinear PDE using Comsol. Here one can find the results u=u(x,y) in the form of an Excel table with the terms {x,y,u}. The domain in which the equation is simulated has the following form enter image description here

Below you can see the solution returned by comsol:

enter image description here

Now, I save this same solution as the Excel file (the reference to it is given above) bring it to Mathematica and build a ListContourPlot:

 ListContourPlot[\[Rho]SNP1, ColorFunction -> "Rainbow", AspectRatio->18/222, ClippingStyle -> Automatic,PlotRange -> {300, 3000}, ImageSize -> 500]

Here is what I see:

enter image description here

One can notice that the geometry of the image I get in Mma differs essentially from the expected one. Namely the narrow stripe part disappeared, and a trapezoidal shape is visible instead. If one inspects the image more closely, he can notice that the yellow part in the Mma image corresponds to the narrow stripe part (that is shown by the red part in the Comsol image). There is, however, a parasitic part of the image, which should not be there at all.

Any idea of how to get rid of it?

$\endgroup$
5
  • 1
    $\begingroup$ ListContourPlot starts by building a Delaunay triangulation of the points, in effect constructing their convex hull. Everything in the convex hull will be included in the plot. We need to find a way not to re-build the triangulation, but instead export the existing one from Comsol, and do the plot based on that. (I know that this comment won't solve the problem, but it's perhaps a start...) $\endgroup$ – Szabolcs Jan 29 '16 at 11:36
  • $\begingroup$ The first question is still: is it possible to export the actual mesh from Comsol (including edges too, rather than only the points)? $\endgroup$ – Szabolcs Jan 29 '16 at 11:41
  • $\begingroup$ @Szabolcs I tried, but did not manage. It is either possible in some special binary format, or as a text file. The latter I did not succeed to transform into expression form. $\endgroup$ – Alexei Boulbitch Jan 29 '16 at 13:39
  • $\begingroup$ Hi Alexei, can you export the *.mphtxt from comsol? Or in in any other way export the connectivity of the mesh elements? $\endgroup$ – user21 Jan 29 '16 at 17:59
  • $\begingroup$ @user21 Hi, I placed the file mesh.mphtxt on the web here: dropbox.com/s/apxvpbiroljetpz/mesh.mphtxt?dl=0 and the file mesh.mphbin here: dropbox.com/s/h257l9t123rieso/mesh.mphbin?dl=0 $\endgroup$ – Alexei Boulbitch Feb 1 '16 at 9:05
4
$\begingroup$

Here is a partial solution. Partial, because it seems to me that the mesh and the solution do not match.

The mphtxt file is a structured representation of a mesh. This file is made up of several sections and comments which start with a # token. There is a coordinates section and then several element sections. To find the starting position of these element sections in the text I look for string tokens in the # section that mark the start in the comments. If I have the position where an element section starts I can then import that section since all element section have the same structure.

Each element section also has several pieces of information. The incidents (ids how the elements are connected to the nodes), markers (each element may have a node attributed). There are also some information I do not know the purpose of, so I just jump over them when reading the file.

A bit more practical: I import the Comsol mphtxt file, delete empty lines and find the position where the coordinate section starts. The start of that section is indicated the # number of mesh points string.

data = DeleteCases[Import["~/Downloads/mesh.mphtxt", "Table"], {}];
startPos = 
 First[Flatten[
   Position[data, {_, "#", "number", "of", "mesh", "points"}]]]
(*15*)

The number before that # tells me how many coordinates there are

data[[startPos]]
(* {282, "#", "number", "of", "mesh", "points"} *)

In this example we have 282 coordinates in the mesh.

nrOfCoords = data[[startPos, 1]];

The next line has information about a convention. In Comsol files it seams that the first node is counted from 0. In Mathematica positions of lists start from 1. So we add an indexOffset of 1 if the file tells us to do so.

temp = data[[startPos + 1, 1]];
indexOffset = If[temp == 0, 1, 0];

Next we read the data from the first coord position until the number of coords that are in the file:

start1 = startPos + 3;
end1 = (startPos + 3) + nrOfCoords - 1;
coords = Developer`ToPackedArray[N[data[[start1 ;; end1]]]];
Developer`PackedArrayQ[coords]
Length[coords] === nrOfCoords
dim = Last[Dimensions[coords]]
(* True *)
(* True *)
(* 2 *)

Now, the element section has a keyword for different element type but besides that each element section has the same structure.

Find the postions of the point elements ("vtx" section), the line elements ("edg" section) and the triangle elements ("tri" section)

pointPos = Flatten[Position[data, {_, "vtx", "#", "type", "name"}]]
(* {302} *)
linePos = Flatten[Position[data, {_, "edg", "#", "type", "name"}]]
(* {326} *)
triPos = Flatten[Position[data, {_, "tri", "#", "type", "name"}]]
(* {766} *)

Since Comsol supports more than the above mentioned elements other files may have different keywords which are not present in this single file. But adding them would not be hard since the structure of those element sections will be the same as the ones which are in this file.

quadPos = Flatten[Position[data, {_, "XXX", "#", "type", "name"}]]
tetPos = Flatten[Position[data, {_, "XXX", "#", "type", "name"}]]
hexPos = Flatten[Position[data, {_, "XXX", "#", "type", "name"}]]

But this is not important to read this specific file.

Next, I wrote a little function to extract the element data. This works again by going through the file finding how many entries need to be read and then read the entries of interest. The minimum we need to get from the element section are the incidents. I also added code to extract the markers. The other entries in the sections I am not sure what they are for, but they are not needed to make a mesh.

getElements[data_, start_, indexOffset_, eleType_] := 
 Block[{pos, nodesPerEle, nrEle, rawIncidents, nrParameters, 
   nrMarkers, markers},
  pos = start + 1;
  (* nodesPerEle=data[[pos,1]]; *)

  (* incidents *)
  pos = pos + 1;
  nrEle = data[[pos, 1]];
  (* jump over the  `# Elements` line *)
  pos = pos + 2;
  rawIncidents = data[[pos ;; pos + nrEle - 1]] + indexOffset;

  (* not sure what parameters are *)
  pos = pos + nrEle + 1;
  nrParameters = data[[pos, 1]];

  (* markers *)
  pos = pos + nrParameters + 2;
  nrMarkers = data[[pos, 1]];
  (* not sure we need/want the marker offset *)
  markers = data[[pos + 2 ;; pos + nrMarkers + 1]] + 1;
  eleType[rawIncidents, Flatten[markers]]
  ]

Now, we actually use that function and get the elements and construct the mesh:

Needs["NDSolve`FEM`"]
pEle = getElements[data, #, indexOffset, PointElement] & /@ pointPos
(*{PointElement[{{2}, {6}, {218}, {238}, {276}, {282}}, {1, 2, 3, 4, 5, 
   6}]}*)

lineEle = 
  getElements[data, #, indexOffset, LineElement] & /@ linePos;
triEle = getElements[data, #, indexOffset, TriangleElement] & /@ 
   triPos;
quadEle = 
 getElements[data, #, indexOffset, TriangleElement] & /@ quadPos
(*{}*)

mesh = ToElementMesh["Coordinates" -> coords, 
  "MeshElements" -> Join[triEle, quadEle], 
  "BoundaryElements" -> lineEle, "PointElements" -> pEle]

There are no quad elements in this mesh. This is just to show what needs to be done in case there were any quad elements. (Add the # token for quad above such that the quadPos would find the start point of that section)

Visualize the wireframe:

mesh["Wireframe"]

enter image description here

We can also have a look at the boundary markers:

mesh["Wireframe"["MeshElement" -> "BoundaryElements", 
  "MeshElementMarkerStyle" -> Blue]]

enter image description here

So far so good. If I import the values at the nodes:

vals = Import["StripeWithNarrowPart1.csv"][[All, 3]];
Length[vals]
(*20938*)

which does not match the number of coordinates in the mesh:

Length[coords]
(*282*)

Either Comsol does not export the true computation mesh or the solution does not match to the mesh. I am sure there is a setting in Comsol to make the mesh and the solution be in sync.

If we can figure this out you should be able to create an interpolating function:

efun = ElementInterpolation[{mesh}, vals]

Hope this helps.

$\endgroup$
7
  • $\begingroup$ @Alexei, did you manage to get the mesh and the values in sync? Also, if you wanted to we could expand this answer a bit by adding more examples and make this an extension of Import. Only if you have the time for this. $\endgroup$ – user21 Feb 2 '16 at 19:00
  • $\begingroup$ I must admit that I do not understand several important things in your answer. I see that it works, but do not understand, how does it work. Probably, partially since I do not understand the content of the mesh file, and the keywords. Or, for example, why do we need to start from 1, rather than from 0. All this prevents me from using it, sorry. $\endgroup$ – Alexei Boulbitch Feb 4 '16 at 10:25
  • $\begingroup$ I added more explanations - hope this helps. Let me know what you do not understand and then I can clarify it. $\endgroup$ – user21 Feb 4 '16 at 18:40
  • $\begingroup$ If anyone else is using this, please @ me. I'd like to extend this but I lack examples (and a comsol license) $\endgroup$ – user21 Feb 14 '16 at 22:08
  • $\begingroup$ @user21, could you have a look at this question... do ou think that this can be done in Mathematica ?stackoverflow.com/questions/43550076/… $\endgroup$ – james Apr 22 '17 at 7:20
4
$\begingroup$

Basing on the advice of @Szabolcs to whom I am very grateful, I think I can offer a workaround.

First I create in Mma the domain in which the equation has been solved in Comsol:

dom = DiscretizeGraphics[Graphics[Polygon[{{0, 0}, {0, 0.018}, {0.155, 
      0.018}, {0.155 + 0.015, 0.008}, {0.222, 0.008}, {0.222, 0}, {0,0}}]]]

looking as follows:

enter image description here

Then basing on the above Excel file I create the interpolation function

  f = Interpolation[\[Rho]SNP1 /. {x_, y_, z_} -> {{x, y}, z}, 
  InterpolationOrder -> 1];

Now it can be plotted within the domain dom:

ContourPlot[f[x, y], {x, y} \[Element] dom, ColorFunction -> "Rainbow", AspectRatio -> 18/222,  ClippingStyle -> Automatic, PlotRange -> {300, 3000}, 
 ImageSize -> 500]

enter image description here

This is already not very bad.

The problem here is that the InterpolationOrder here can only be 1, since the array returned by Comsol is not structured. The data is, thus, a bit rougher than the solution allows.

May be, you will propose a better way?

$\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.