MEDCoupling : problem using field (with profile) and buildSubPart

Development
1

Hi!
I am trying to read a subpart (defined by a mesh group) of a field, using MEDCoupling.
While following some examples, I could make it work if the field is defined on the whole mesh, but if the field is defined using a profile, buildSubPart throws an error : “medcoupling.InterpKernelException: DataArrayDouble::selectByTupleIdSafe : some ids has been detected to be out of [0,this->getNumberOfTuples) !”

I attached a simple python with a simple MED file.

Here is the python code:

from medcoupling import *

fname="box_profile.med"
meshName="mesh"
fieldName="reslin__DEPL"
iteration=1
order=1

medfileField1TS=MEDFileField1TS.New(fname,fieldName,iteration,order)
mm=MEDFileMesh.New(fname)
valsRead,pflRead=medfileField1TS.getFieldWithProfile(ON_NODES,0,mm)
valsRead.setName("")
pflRead.setName("PROFIL__00000002")

group_cell_ids = mm.getGroupArr(0, "G1", renum=False)
print(f"{group_cell_ids=}")

whole_mesh = mm.getMeshAtLevel(0)
f1 = MEDCouplingFieldDouble.New(ON_NODES, ONE_TIME)
f1.setName(fieldName)
f1.setMesh(whole_mesh)
f1.setArray(valsRead)
f1.setTime(0.0, iteration, order)

f2=f1.buildSubPart(group_cell_ids) # InterpKernelException here

Could you please help me understand what is going wrong here ?
Thank you in advance!
All the best,

test_group_profile.py (656 Bytes)
box_profile.med (38.1 KB)

2

I think I found some way to make it work, I will leave it here in case it will be useful to someone else and hoping that a kind soul can to point me towards a simpler&better way… :slight_smile:

from medcoupling import *

fname="box_profile.med"
meshName="mesh"
fieldName="reslin__DEPL"
iteration=1
order=1

medfileField1TS=MEDFileField1TS.New(fname,fieldName,iteration,order)
mm=MEDFileMesh.New(fname)
valsRead,pflRead=medfileField1TS.getFieldWithProfile(ON_NODES,0,mm)

whole_mesh = mm.getMeshAtLevel(0)
profile_cell_ids = whole_mesh.getCellIdsLyingOnNodes(pflRead, fullyIn=True)

computed_mesh=whole_mesh[profile_cell_ids]
computed_mesh.zipCoords()
computed_mesh.setName(mm.getName())

f1 = MEDCouplingFieldDouble.New(ON_NODES, ONE_TIME)
f1.setName(fieldName)
f1.setMesh(computed_mesh)
f1.setArray(valsRead)

group_cell_ids = mm.getGroupArr(0, "G1", renum=False)
group_cellids_in_profile = group_cell_ids.buildIntersection(profile_cell_ids)

group_mesh = whole_mesh[group_cellids_in_profile]

node_ids_new, num_new_node_ids = group_mesh.getNodeIdsInUse()
pflRead_group = node_ids_new.invertArrayO2N2N2O(num_new_node_ids)
pflRead_group.setName("PROFIL__00000002")
f2 = f1.buildSubPart(group_cellids_in_profile)


ff=MEDFileField1TS.New()
ff.setFieldProfile(f2,mm,0,pflRead_group)

outfname = "box_profile_group.med"
mm.write(outfname,2)
ff.write(outfname,0)

box_profile_group.med (35.0 KB)
test_group_profile.py (1.1 KB)

3

Hi Luca,

You got it right.

G1 is made of cells [3, 4], whereas the profile on nodes gives the cells [0, 1, 4, 5, 7]. 3 is missing, which gives you the initial error.

Doing an intersection of the two dataarrays like you did is the good solution. Good job :+1:

All the best,

Christophe

4

Hi Christophe,

Thank you for looking into this and giving me your feedback! It feels much better to know that I am on the right way. In case someone is interested, I am writing a medcoupling wrapper to make classic civil engineering post-processing easier (I need it for myself), everything is opensourced. It is still at a very basic stage, if you want to have a look at it I put it here:

If you have any idea or advice, I am all ears ! :slight_smile:

Have a nice day, all the best,
Luca

2 Likes