[Gmsh] Restrict anisotropic mesh

[Lionel Wilhelm]

Dear all,

I'm trying to restrict an anisotropic field to a given region but the 
command Restrict does not seem be compliant with anisotropic fields. The 
restricted field has the MAX_LC value everywhere. I had a look in the 
source code 
(https://gitlab.onelab.info/gmsh/gmsh/blob/903725dc22de1f2e6545d4f843588f74df5492f4/Mesh/Field.cpp, 
line 1771) but I do not see the problem as the RestrictField class only 
deals with entity id's and fields (whether isotropic or anisotropic)

In order to illustrate, let me join a sample .geo file. The 
FieldToRestrict variable allow to restrict either an isotropic or an 
anisotropic mesh.

Which workaround would you use ? would that be hard to implement ? would 
you have some advices ?

Kindly,

Lionel Wilhelm

-------------- next part --------------
SetFactory("OpenCASCADE");

Mesh.Algorithm = 7; //Bamg
Mesh.Algorithm3D= 7; //mmg3d

// Create an empty list to store the individual boxes id.
boxes = {};

//create a new box
boxes += newv;
Box(boxes[0]) = {0,0,0,1,1,1};

// anisotropic field definition
Field_number = 1;
Field[Field_number] = MathEvalAniso;   
Field[Field_number].m11 = "1/(0.1)^2";
Field[Field_number].m12 = "0";
Field[Field_number].m13 = "0";
Field[Field_number].m22 = "1/(0.1)^2";
Field[Field_number].m23 = "0";
Field[Field_number].m33 = "1/(0.01)^2";

// Isotropic field to restrict
Field_number = 2;
Field[Field_number] = MathEval;
Field[Field_number].F = "0.05";

// Restriction
FieldToRestrict = 1; // 1=anisotropic field, 2 = isotropic field

region_id = boxes[0];
region_faces = Unique(Abs(Boundary{Volume{region_id[]};}));
region_edges = Unique(Abs(Boundary{Surface{region_faces[]};}));
region_points= Unique(Abs(Boundary{Line{region_edges[]};}));
Field_number = 3;
Field[Field_number] = Restrict;
Field[Field_number].IField = FieldToRestrict;
Field[Field_number].RegionsList = {region_id[]};
Field[Field_number].VerticesList = {region_points[]};
Field[Field_number].EdgesList = {region_edges[]};
Field[Field_number].FacesList = {region_faces[]};

Background Field = 3;