Unexpected view of the Comsol result when I postprocess it in Mma

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:

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]

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?

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"]

We can also have a look at the boundary markers:

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

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.