I am trying to do a finiteElement Simulation of the Sun - Earth system for some modified gravity model in cylindrical coordinates. I find that 2nd order elements approximate the gravitational field better close to the Earth and first order Elements better close to the sun and would therefore like to have a mesh with varying element order. Is something like this possible?

An example code of roughly what I would like to get is here (however unforuntaley using first order elements everywhere):



cutoff = 400; (*Radius of the sun is 1 in the used units*)

A = ImplicitRegion[True, {{x, 0, cutoff}, {y, 0, cutoff}}];

f = Function[{vertices, 
    area}, (area > 
      0.001*(1 + 
         10*Abs[((Mean[vertices][[1]])^2 + (Mean[vertices][[2]])^2 - 
    (area > 
       0.00005 && (0.8^2 <= (Mean[vertices][[
           1]])^2 + (Mean[vertices][[2]])^2 <= 1.1))
    (area > 
      0.0005^2*(1 + 
         10^4*Abs[((Mean[vertices][[1]])^2 + (Mean[vertices][[2]] - 
               dist)^2 - REN^2)]))
    || (area > 
       0.05 && (0 <= (Mean[vertices][[1]])^2 + (Mean[vertices][[
           2]])^2 <= (0.8 REN)^2))];

B = DiscretizeRegion[A, MeshRefinementFunction -> f];

domb = MeshOrderAlteration[B["MakeRepresentation"["ElementMesh"]], 1];

I tried to create a custom mesh by manually merging two different meshes roughly like this (splitting the domain in two and using first order elements for one part, second order for the other part, points contains the points of both domains, Elements and Elements2 the information how they should be combined to elements):

dombFinal = 
  ToElementMesh["Coordinates" -> points, 
   "MeshElements" -> {TriangleElement[Elements], 

Where Elements are 1st order Elements and Elements2 second order Elements but I get error messages that my element order is not consistent. Is there a workaround? Thanks in advance!


1 Answer 1


No, you can not have mixed element order in one and the same a mesh. The only thing you can have, is mixed order for multiple dependent variables. Say, your PDE consists of u, v and p, then you can use InterpolationOrder to specify that, for example, u and v are second order, while p is first order. But you can not have u first order in some spatial part of the mesh and second order in another.

  • $\begingroup$ Thank you for the answer. That is unforunate but since it appears to be the correct answer I marked is as such. $\endgroup$
    – Eddi
    Feb 15, 2022 at 9:23
  • $\begingroup$ @Eddi, thanks. Have you ever seen a FEM system that actually allows for something like what you are asking for? If so, could you send me a link to that. Currently, I think, it's technically not possible to do this but I might be wrong. I think, adaptive mesh refinement, is more what you'd use in such a scenario. (but V13.0 does not have that either) Another idea, maybe just maybe, would be to use ToGradedMesh. $\endgroup$
    – user21
    Feb 15, 2022 at 9:37
  • $\begingroup$ Yes something like this is possible indeed and goes by the name hp-FEM, en.wikipedia.org/wiki/Hp-FEM. I learned this in my FEM class and I know that the software NGSolve (ngsolve.org) is capable of this. By default FEM spaces are defined by the maximum order there and contain all polynomials of equal or smaller order. This allows for something called hp-refinement, where you can locally make a triangle smaller (h refinement) or increase the polynomial degree (p refinement). The documentation is not that great unfortunately. $\endgroup$
    – Eddi
    Feb 16, 2022 at 6:53
  • $\begingroup$ Thanks for the link - I did not get that you are looking for hp-adaptivity. $\endgroup$
    – user21
    Feb 16, 2022 at 19:02

Your Answer

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

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