# Finite Element Mesh generation with varying element order

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):

Needs["NDSolveFEM"];

dist=215.0321556705476;
REN=0.009167888457668535;

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][])^2 + (Mean[vertices][])^2 -
1^2)]))||
(area >
0.00005 && (0.8^2 <= (Mean[vertices][[
1]])^2 + (Mean[vertices][])^2 <= 1.1))
||
(area >
0.0005^2*(1 +
10^4*Abs[((Mean[vertices][])^2 + (Mean[vertices][] -
dist)^2 - REN^2)]))
|| (area >
0.05 && (0 <= (Mean[vertices][])^2 + (Mean[vertices][[
2]])^2 <= (0.8 REN)^2))];

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

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


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],
TriangleElement[Elements2]}];


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!

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.