MEDCoupling C++ examples

Fields

Create

Standard build of a tensor field on cells with no time attached

  MEDCoupling::MEDCouplingFieldDouble* fieldOnCells=MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_CELLS,MEDCoupling::NO_TIME);
  fieldOnCells->setName("MyTensorFieldOnCellNoTime");
  fieldOnCells->setMesh(mesh);
  mesh->decrRef(); // no more need of mesh because mesh has been attached to fieldOnCells
  MEDCoupling::DataArrayDouble *array=MEDCoupling::DataArrayDouble::New();
  array->alloc(fieldOnCells->getMesh()->getNumberOfCells(),9);//Implicitly fieldOnCells will be a 9 components field.
  array->fillWithValue(7.);
  fieldOnCells->setArray(array);
  array->decrRef();
  // fieldOnCells is now usable
  // ...
  // fieldOnCells is no more useful here : release it
  fieldOnCells->decrRef();

Standard build of a scalar field on nodes with no time attached

  MEDCoupling::MEDCouplingFieldDouble* fieldOnNodes=MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_NODES,MEDCoupling::NO_TIME);
  fieldOnNodes->setName("MyScalarFieldOnNodeNoTime");
  fieldOnNodes->setMesh(mesh);
  mesh->decrRef(); // no more need of mesh because mesh has been attached to fieldOnNodes
  MEDCoupling::DataArrayDouble *array=MEDCoupling::DataArrayDouble::New();
  array->alloc(fieldOnNodes->getMesh()->getNumberOfNodes(),1);//Implicitly fieldOnNodes will be a 1 component field.
  array->fillWithValue(8.);
  fieldOnNodes->setArray(array);
  array->decrRef();
  // fieldOnNodes is now usable
  // ...
  // fieldOnNodes is no more useful here : release it
  fieldOnNodes->decrRef();

Standard build of a vector field on cells with one time attached and no time interval

  MEDCoupling::MEDCouplingFieldDouble* fieldOnCells=MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_CELLS,MEDCoupling::ONE_TIME);
  fieldOnCells->setName("MyTensorFieldOnCellNoTime");
  fieldOnCells->setTimeUnit("ms"); // Time unit is ms.
  fieldOnCells->setTime(4.22,2,-1); // Time attached is 4.22 ms, iteration id is 2 and order id (or sub iteration id) is -1
  fieldOnCells->setMesh(mesh);
  mesh->decrRef(); // no more need of mesh because mesh has been attached to fieldOnCells
  MEDCoupling::DataArrayDouble *array=MEDCoupling::DataArrayDouble::New();
  array->alloc(fieldOnCells->getMesh()->getNumberOfCells(),2);//Implicitly fieldOnCells will be a 2 components field.
  array->fillWithValue(7.);
  fieldOnCells->setArray(array);
  array->decrRef();
  // fieldOnCells is now usable
  // ...
  // fieldOnCells is no more useful here : release it
  fieldOnCells->decrRef();

Standard build of a vector field on nodes defined on a time interval with a constant value during this interval

  MEDCoupling::MEDCouplingFieldDouble* fieldOnNodes=MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_NODES,MEDCoupling::CONST_ON_TIME_INTERVAL);
  fieldOnNodes->setName("MyVecFieldOnNodeWithConstTime");
  fieldOnNodes->setTimeUnit("ms"); // Time unit is ms.
  fieldOnNodes->setStartTime(4.22,2,-1);
  fieldOnNodes->setEndTime(6.44,4,-1); // fieldOnNodes is defined in interval [4.22 ms,6.44 ms] 
  fieldOnNodes->setMesh(mesh);
  mesh->decrRef(); // no more need of mesh because mesh has been attached to fieldOnNodes
  MEDCoupling::DataArrayDouble *array=MEDCoupling::DataArrayDouble::New();
  array->alloc(fieldOnNodes->getMesh()->getNumberOfNodes(),3);//Implicitly fieldOnNodes will be a 3 components field.
  array->fillWithValue(8.);
  fieldOnNodes->setArray(array);
  array->decrRef();
  // fieldOnNodes is now usable
  // ...
  // fieldOnNodes is no more useful here : release it
  fieldOnNodes->decrRef();

Getting maximal and minimal fields

In this example we

• create two fields with two tuples per two components,
• use MaxFields() to get a field holding maximal values of the two fields.
• use MinFields() to get a field holding minimal values of the two fields.

  const double vals1[4]   = {0.,2., 4.,6.}; // for field 1
  const double vals2[4]   = {2.,0., 6.,4.}; // for field 2
  const double valsMax[4] = {2.,2., 6.,6.}; // expected max field
  const double valsMin[4] = {0.,0., 4.,4.}; // expected min field
  // field 1
  MCAuto<DataArrayDouble> valsArr1 = DataArrayDouble::New();
  valsArr1->useExternalArrayWithRWAccess( vals1, 2,2 ); // 2 tuples per 2 components
  MCAuto<MEDCouplingFieldDouble> field1 = MEDCouplingFieldDouble::New( ON_NODES );
  field1->setArray( valsArr1 );
  // field 2
  MCAuto<DataArrayDouble> valsArr2 = DataArrayDouble::New();
  valsArr2->useExternalArrayWithRWAccess( vals2, 2,2 ); // 2 tuples per 2 components
  MCAuto<MEDCouplingFieldDouble> field2 = MEDCouplingFieldDouble::New( ON_NODES );
  field2->setArray( valsArr2 );
  // max field 
  MCAuto<MEDCouplingFieldDouble> fieldMax = MEDCouplingFieldDouble::MaxFields( field1, field2 );
  CPPUNIT_ASSERT( std::equal( valsMax, valsMax+4, fieldMax->getArray()->getConstPointer() )); // fieldMax == valsMax
  // min field 
  MCAuto<MEDCouplingFieldDouble> fieldMin = MEDCouplingFieldDouble::MinFields( field1, field2 );
  CPPUNIT_ASSERT( std::equal( valsMin, valsMin+4, fieldMin->getArray()->getConstPointer() )); // fieldMin == valsMin

Concatenating fields

In this example we

• create an 1D mesh and a field on it,
• make a deep copy of the mesh and the field,
• translate the mesh and the field,
• use two variants of MergeFields() to create one field from the two by concatenating them and their meshes.

  // mesh1
  const double coords[3] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh1 = MEDCouplingCMesh::New();
  mesh1->setCoords(coordsArr); // mesh becomes a 1D
  // field1
  MCAuto<MEDCouplingFieldDouble> field1 =
    mesh1->fillFromAnalytic( ON_CELLS, 1, "x");
 
  // mesh2 and field2
  MCAuto<MEDCouplingFieldDouble> field2 =
    field1->cloneWithMesh( true );
  double vec[1] = { 5. };
  (const_cast<MEDCoupling::MEDCouplingMesh *>(field2->getMesh()))->translate(vec); // translate mesh2
  field2->applyFunc("x + 5"); // "translate" field2
 
  // concatenate field1 and field2
  MCAuto<MEDCouplingFieldDouble> field3 =
    MEDCouplingFieldDouble::MergeFields( field1, field2 );
  std::vector<const MEDCouplingFieldDouble *> fields( 2 );
  fields[0] = field1;
  fields[1] = field2;
  MCAuto<MEDCouplingFieldDouble> field4 =
    MEDCouplingFieldDouble::MergeFields( fields );

The result field is twice "longer" than field1.

Getting a field copy with different time discretization

First, we create a supporting 2D mesh and a field on it got using fillFromAnalytic(). Time discretization of this field is ONE_TIME.

  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh = MEDCouplingCMesh::New();
  mesh->setCoords(coordsArr,coordsArr);
  MCAuto<MEDCouplingFieldDouble> field1 =
    mesh->fillFromAnalytic( MEDCoupling::ON_NODES,1,"x+y");
  CPPUNIT_ASSERT( field1->getTimeDiscretization() == MEDCoupling::ONE_TIME );

Now we use buildNewTimeReprFromThis() to get a copy of field1 whose time discretization is NO_TIME.

  MCAuto<MEDCouplingFieldDouble> field2 =
    field1->buildNewTimeReprFromThis( MEDCoupling::NO_TIME, false );
  CPPUNIT_ASSERT( field2->getTimeDiscretization() == MEDCoupling::NO_TIME );

Creating a field using an expression

First, we create a 2D Cartesian mesh constituted by 2 cells.

  const double coords[4] = {0.,2.,4.,6.}; // 6. is not used
  MCAuto<DataArrayDouble> x = DataArrayDouble::New();
  x->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<DataArrayDouble> y = DataArrayDouble::New();
  y->useExternalArrayWithRWAccess( coords, 2, 1 );
  MCAuto<MEDCouplingCMesh> mesh=MEDCouplingCMesh::New();
  mesh->setCoords(x,y);

Now we use fillFromAnalytic3() to get a MEDCouplingFieldDouble on cells filled with values computed using an expression func. This expression is applied to coordinates of each point (barycenter) for which the field value is computed. We want to get the field on cells, with 3 components computed as follows. (In func, we refer to the first component of a point using the variable "a", and to the second component, using the variable "b").

• Component #0 = the second coordinate of the point hence "IVec * b" in func.
• Component #1 = the first coordinate of the point hence "JVec * a".
• Component #2 = distance between the point and SC origin (0.,0.) hence "KVec * sqrt( a*a + b*b )".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  const char func[] = "IVec * b + JVec * a + KVec * sqrt( a*a + b*b ) + 10";
  const char* varNames[2] = { "a", "b" }; // names used to refer to X and Y coord components
  std::vector<std::string> varNamesVec( varNames, varNames+2 );
  MCAuto<MEDCouplingFieldDouble> field =
    mesh->fillFromAnalyticNamedCompo( MEDCoupling::ON_CELLS, 3, varNamesVec, func );

Now we ascertain that the result field is as we expect. We check the second tuple of the field. We get barycenter of the cell #1 and checks that values of the second tuple are computed as we want.

  double val1[3]; // a value (vector) of the cell #1
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 ); // 3 components in the field
  field->getArray()->getTuple( 1, val1 );
  //
  MCAuto<DataArrayDouble> bc =
    mesh->computeCellCenterOfMass(); // func is applied to barycenters of cells
  double bc1[2]; // coordinates of the second point
  bc->getTuple( 1, bc1 );
  //
  double dist = sqrt( bc1[0]*bc1[0] + bc1[1]*bc1[1] );  // "sqrt( a*a + b*b )"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[0], 10 + bc1[1], 13 ); // "10 + IVec * b"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[1], 10 + bc1[0], 13 ); // "10 + JVec * a"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[2], 10 + dist  , 13 ); // "10 + KVec * sqrt( a*a + b*b )"

Creating a field using an expression

First, we create a 2D Cartesian mesh constituted by 2 cells. Note that we set names to coordinates arrays ("a" and "b" ) which will be used to refer to corresponding coordinates within a function.

  const double coords[4] = {0.,2.,4.,6.}; // 6. is not used
  MCAuto<DataArrayDouble> x = DataArrayDouble::New();
  x->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<DataArrayDouble> y = DataArrayDouble::New();
  y->useExternalArrayWithRWAccess( coords, 2, 1 );
  x->setInfoOnComponent(0,"a"); //  name used to refer to X coordinate within a function
  y->setInfoOnComponent(0,"b"); //  name used to refer to Y coordinate within a function
  MCAuto<MEDCouplingCMesh> mesh=MEDCouplingCMesh::New();
  mesh->setCoords(x,y);

Now we use fillFromAnalytic2() to get a MEDCouplingFieldDouble on cells filled with values computed using an expression func. This expression is applied to coordinates of each point (barycenter) for which the field value is computed. We want to get the field on cells, with 3 components computed as follows. (In func, we refer to the first component of a point using the variable "a", and to the second component, using the variable "b").

• Component #0 = the second coordinate of the point hence "IVec * b" in func.
• Component #1 = the first coordinate of the point hence "JVec * a".
• Component #2 = distance between the point and SC origin (0.,0.) hence "KVec * sqrt( a*a + b*b )".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  const char func[] = "IVec * b + JVec * a + KVec * sqrt( a*a + b*b ) + 10";
  MCAuto<MEDCouplingFieldDouble> field =
    mesh->fillFromAnalyticCompo( MEDCoupling::ON_CELLS, 3, func );

Now we ascertain that the result field is as we expect. We check the second tuple of the field. We get barycenter of the cell #1 and checks that values of the second tuple are computed as we want.

  double val1[3]; // a value (vector) of the cell #1
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 ); // 3 components in the field
  field->getArray()->getTuple( 1, val1 );
  //
  MCAuto<DataArrayDouble> bc =
    mesh->computeCellCenterOfMass(); // func is applied to barycenters of cells
  double bc1[2]; // coordinates of the second point
  bc->getTuple( 1, bc1 );
  //
  double dist = sqrt( bc1[0]*bc1[0] + bc1[1]*bc1[1] );  // "sqrt( a*a + b*b )"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[0], 10 + bc1[1], 13 ); // "10 + IVec * b"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[1], 10 + bc1[0], 13 ); // "10 + JVec * a"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[2], 10 + dist  , 13 ); // "10 + KVec * sqrt( a*a + b*b )"

Creating a field using an expression

First, we create a 2D Cartesian mesh constituted by 2 cells.

  const double coords[4] = {0.,2.,4.,6.}; // 6. is not used
  MCAuto<DataArrayDouble> x = DataArrayDouble::New();
  x->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<DataArrayDouble> y = DataArrayDouble::New();
  y->useExternalArrayWithRWAccess( coords, 2, 1 );
  MCAuto<MEDCouplingCMesh> mesh=MEDCouplingCMesh::New();
  mesh->setCoords(x,y);

Now we use fillFromAnalytic() to get a MEDCouplingFieldDouble on cells filled with values computed using an expression func. This expression is applied to coordinates of each point (barycenter) for which the field value is computed. We want to get the field on cells, with 3 components computed as follows. (In func, we refer to the first component of a point using the variable "a", and to the second component, using the variable "b").

• Component #0 = the second coordinate of the point hence "IVec * b" in func.
• Component #1 = the first coordinate of the point hence "JVec * a".
• Component #2 = distance between the point and SC origin (0.,0.) hence "KVec * sqrt( a*a + b*b )".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  const char func[] = "IVec * b + JVec * a + KVec * sqrt( a*a + b*b ) + 10";
  MCAuto<MEDCouplingFieldDouble> field =
    mesh->fillFromAnalytic( MEDCoupling::ON_CELLS, 3, func );

Now we ascertain that the result field is as we expect. We check the second tuple of the field. We get barycenter of the cell #1 and checks that values of the second tuple are computed as we want.

  double val1[3]; // a value (vector) of the cell #1
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 ); // 3 components in the field
  field->getArray()->getTuple( 1, val1 );
  //
  MCAuto<DataArrayDouble> bc =
    mesh->computeCellCenterOfMass(); // func is applied to barycenters of cells
  double bc1[2]; // coordinates of the second point
  bc->getTuple( 1, bc1 );
  //
  double dist = sqrt( bc1[0]*bc1[0] + bc1[1]*bc1[1] );  // "sqrt( a*a + b*b )"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[0], 10 + bc1[1], 13 ); // "10 + IVec * b"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[1], 10 + bc1[0], 13 ); // "10 + JVec * a"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[2], 10 + dist  , 13 ); // "10 + KVec * sqrt( a*a + b*b )"

Creation of a sub part of a field

Creation of a sub part of a field on cells

  MEDCoupling::MEDCouplingUMesh *mesh1=MEDCoupling::MEDCouplingBasicsTest::build2DTargetMesh_1();
  MEDCoupling::MEDCouplingFieldDouble *f1=MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_CELLS,MEDCoupling::ONE_TIME);
  f1->setTime(2.3,5,6);
  f1->setMesh(mesh1);
  MEDCoupling::DataArrayDouble *array=MEDCoupling::DataArrayDouble::New();
  array->alloc(mesh1->getNumberOfCells(),2);
  const double arr1[10]={3.,103.,4.,104.,5.,105.,6.,106.,7.,107.};
  std::copy(arr1,arr1+10,array->getPointer());
  f1->setArray(array);
  array->decrRef();

The field on cells f1 lies on a mesh containing 5 cells and 9 nodes. So this field f1 contains 5 tuples of 2 components each (10 values). Now let's create a subfield on cells f2 from f1.

  const mcIdType part1[3]={2,1,4};
  MEDCoupling::MEDCouplingFieldDouble *f2=f1->buildSubPart(part1,part1+3);

f1 is a field on cells, buildSubPart method performs an extraction on cells too.

So the array part1 lists ids on cells.

• cell #0 of f2 is the same cell of cell #2 of f1
• cell #1 of f2 is the same cell of cell #1 of f1
• cell #2 of f2 is the same cell of cell #4 of f1

So f2 contains 3 tuples with 2 components.

The underlying mesh of f2 contains a newly created mesh with 3 cells (not as mesh1 in f1) and 9 nodes (as mesh1 in f1).
For fields on cells the number of tuples of the returned field is always equal to the number of ids given in input (here part1).
Only fields on cells have this particular behaviour.

Creation of a sub part of a field on nodes

  f1=MEDCoupling::MEDCouplingFieldDouble::New(MEDCoupling::ON_NODES,MEDCoupling::ONE_TIME);
  f1->setTime(2.3,5,6);
  f1->setMesh(mesh1);
  array=MEDCoupling::DataArrayDouble::New();
  array->alloc(mesh1->getNumberOfNodes(),2);
  const double arr2[18]={3.,103.,4.,104.,5.,105.,6.,106.,7.,107.,8.,108.,9.,109.,10.,110.,11.,111.};
  std::copy(arr2,arr2+18,array->getPointer());  
  f1->setArray(array);
  array->decrRef();

The field on nodes f1 lies on a mesh containing 5 cells and 9 nodes. So this field f1 contains 9 tuples of 2 components each (18 values). Now let's create a subfield on nodes f2 from f1.

  const mcIdType part2[2]={1,2};
  f2=f1->buildSubPart(part2,part2+2);

f1 is a field on nodes, but buildSubPart method performs an extraction on cells.

After the call of buildSubPart on node field f1, f1 will be reduced on a submesh of mesh1 containing cells whose ids are in part2. So here the number of cells of f2 is 2 and the number of nodes is 4.
So contrary to fields on cells, it is normal for fields on nodes that number of tuples of the returned field of buildSubPart method does not match the size of the input array (here part2).

Modify

Subtracting field on different meshes

We make two meshes in 1D space with no cells and 4 nodes. Nodes #0 and #2 are swapped in the two meshes.
And we make two fields on these meshes, so that fields values to equal to node coordinates of the underlying meshes.

  const double coords1[4] = {0.,1.,2.,3.};
  const double coords2[4] = {2.,1.,0.,3.}; //  #0 <==> #2
  // mesh 1
  MCAuto<MEDCouplingUMesh> mesh1 = MEDCouplingUMesh::New();
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords1, 4, 1 );
  mesh1->setCoords(coordsArr);
  mesh1->setMeshDimension(0);
  mesh1->allocateCells(0);
  mesh1->finishInsertingCells();
  // mesh 2
  MCAuto<MEDCouplingUMesh> mesh2 =
    (MEDCouplingUMesh*) mesh1->deepCopy();
  mesh2->getCoords()->useExternalArrayWithRWAccess( coords2, 4, 1 );

We are going to subtract field2 from field1, though they are on different meshes. substractInPlaceDM() allows us doing this. We use a mesh comparison level levOfCheck = 10 that allows subtracting fields on meshes with different node arrays.

  MCAuto<MEDCouplingFieldDouble> field1 =
    mesh1->fillFromAnalytic( MEDCoupling::ON_NODES,1,"x"); // field1 values == coords1
  MCAuto<MEDCouplingFieldDouble> field2 =
    mesh2->fillFromAnalytic( MEDCoupling::ON_NODES,1,"x"); // field2 values == coords2
  const int levOfCheck = 10; // nodes can be permuted
  field1->substractInPlaceDM( field2, levOfCheck, 1e-13, 0 ); // values #0 and #2 must swap

After applying substractInPlaceDM() the both fields lie on mesh2. As substractInPlaceDM() permutes values of field1 before value subtraction, and thus field1 becomes equal to field2, hence their subtraction results in a zero field.

  field2->applyFunc( 1, 0.0 ); // all field2 values == 0.0
  CPPUNIT_ASSERT( field1->isEqual( field2, 1e-13, 1e-13 )); // field1 == field2 == 0.0

Changing the underlying mesh

We make two meshes in 1D space with no cells and 4 nodes. Nodes #0 and #2 are swapped in the two meshes.

  const double coords1[4] = {0.,1.,2.,3.};
  const double coords2[4] = {2.,1.,0.,3.}; //  #0 <==> #2
  // mesh 1
  MCAuto<MEDCouplingUMesh> mesh1 = MEDCouplingUMesh::New();
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords1, 4, 1 );
  mesh1->setCoords(coordsArr);
  mesh1->setMeshDimension(0);
  mesh1->allocateCells(0);
  mesh1->finishInsertingCells();
  // mesh 2
  MCAuto<MEDCouplingUMesh> mesh2 =
    (MEDCouplingUMesh*) mesh1->deepCopy();
  mesh2->getCoords()->useExternalArrayWithRWAccess( coords2, 4, 1 );

We are going to use changeUnderlyingMesh() to set mesh2 instead of mesh1 as a support of a field.
We use fillFromAnalytic() to make a field on nodes of mesh1, so that its values to equal to node coordinates. Then we use changeUnderlyingMesh() to change the underlying mesh of the field. (We use a mesh comparison level levOfCheck = 10 that allows substituting meshes with different node arrays.) As a result, we expect that values of the field are also permuted same as nodes of the two meshes, and thus its values become equal to the array coords2.

  MCAuto<MEDCouplingFieldDouble> field =
    mesh1->fillFromAnalytic( MEDCoupling::ON_NODES,1,"x"); // field values == coords1
  const int levOfCheck = 10; // nodes can be permuted
  field->changeUnderlyingMesh( mesh2, levOfCheck, 1e-13, 0 ); // values #0 and #2 must swap
  CPPUNIT_ASSERT( std::equal( coords2, coords2+4, field->getArray()->getConstPointer() ));

Changing a field using an expression

We create a 2D vector field with 2 tuples and we want to transform this field using an expression using applyFunc(). The expression func is applied to each atomic value of the field. We want to change the field as follows. (In func, we use the variable "v" to refer to an atomic field value).

• Component #0 = component #0 (remains the same) hence "IVec * v" in func.
• Component #1 = component #1 ^ 2 hence "JVec * v*v".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  const double v[4] = {1.,2., 3.,4.};
  MCAuto<DataArrayDouble> array = DataArrayDouble::New();
  array->useExternalArrayWithRWAccess( v, 2, 2 ); // 2 tuples per 2 components
  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS );
  field->setArray( array );
  const char func[] = "IVec * v + JVec * w*w + 10";
  field->applyFunc( 2, func );
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 2 ); // 2 components remains

Now we ascertain that the result field is as we expect.

  const double* v2 = field->getArray()->getConstPointer();
  CPPUNIT_ASSERT_DOUBLES_EQUAL( v2[0], 10 + v[0], 13 );      // "10 + IVec * v"  
  CPPUNIT_ASSERT_DOUBLES_EQUAL( v2[1], 10 + v[1]*v[1], 13 ); // "10 + JVec * v*v"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( v2[2], 10 + v[2], 13 );      // "10 + IVec * v"  
  CPPUNIT_ASSERT_DOUBLES_EQUAL( v2[3], 10 + v[3]*v[3], 13 ); // "10 + JVec * v*v"

Changing a field using an expression

We create a 2D vector field with 2 values (vectors) and then we transform this field into a 3D vector field by applying an expression to values of the 2D field using applyFunc3(). The expression func is applied to components of each vector of the field. We want the field to have 3 components computed as follows. (In func, we refer to the first component of a field value using the variable "a", and to the second component, using the variable "b", as we define it by varNamesVec).

• Component #0 = the second vector component hence "IVec * b" in func.
• Component #1 = the first vector component hence "JVec * a".
• Component #2 = a vector magnitude hence "KVec * sqrt( a*a + b*b )".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  // create a 2D vector field
  const double values[4] = {1.,1., 2.,1.};
  MCAuto<DataArrayDouble> array = DataArrayDouble::New();
  array->useExternalArrayWithRWAccess( values, 2, 2 ); // 2 tuples per 2 components
  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS );
  field->setArray( array );
  // transform the field to a 3D vector field
  const char func[] = "IVec * b + JVec * a + KVec * sqrt( a*a + b*b ) + 10";
  const char* varNames[2] = { "a", "b" }; // names used to refer to X and Y components
  std::vector<std::string> varNamesVec( varNames, varNames+2 );
  field->applyFuncNamedCompo( 3, varNamesVec, func ); // require 3 components 
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 ); // 3 components as required

Now we ascertain that the result field is as we expect. We check the second vector of the field.

  double vec1[3]; // vector #1
  field->getArray()->getTuple( 1, vec1 );
  const double a = values[2], b = values[3]; // initial components of the vector #1
  CPPUNIT_ASSERT_DOUBLES_EQUAL( vec1[0], 10 + b, 13 ); // "10 + IVec * b"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( vec1[1], 10 + a, 13 ); // "10 + JVec * a"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( vec1[2], 10 + sqrt(a*a+b*b), 13 ); // "10 + KVec * sqrt( a*a + b*b )"

Changing a field using an expression

We create a 2D vector field with 2 values (vectors) and then we transform this field into a 3D vector field by applying an expression to values of the 2D field using applyFunc2(). Note that we set component info the array ("a" and "b" ) which will be used to refer to corresponding components within a function. The expression func is applied to components of each vector of the field. We want the field to have 3 components computed as follows. (In func, we refer to the first component of a field value using the variable "a", and to the second component, using the variable "b").

• Component #0 = the second vector component hence "IVec * b" in func.
• Component #1 = the first vector component hence "JVec * a".
• Component #2 = a vector magnitude hence "KVec * sqrt( a*a + b*b )".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  // create a 2D vector field
  const double values[4] = {1.,1., 2.,1.};
  MCAuto<DataArrayDouble> array = DataArrayDouble::New();
  array->useExternalArrayWithRWAccess( values, 2, 2 ); // 2 tuples per 2 components
  array->setInfoOnComponent(0,"a"); // name used to refer to X component within a function
  array->setInfoOnComponent(1,"b"); // name used to refer to Y component within a function
  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS );
  field->setArray( array );
  // transform the field to a 3D vector field
  const char func[] = "IVec * b + JVec * a + KVec * sqrt( a*a + b*b ) + 10";
  field->applyFuncCompo( 3, func ); // require 3 components 
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 ); // 3 components as required

Now we ascertain that the result field is as we expect. We check the second vector of the field.

  double vec1[3]; // vector #1
  field->getArray()->getTuple( 1, vec1 );
  const double a = values[2], b = values[3]; // initial components of the vector #1
  CPPUNIT_ASSERT_DOUBLES_EQUAL( vec1[0], 10 + b, 13 ); // "10 + IVec * b"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( vec1[1], 10 + a, 13 ); // "10 + JVec * a"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( vec1[2], 10 + sqrt(a*a+b*b), 13 ); // "10 + KVec * sqrt( a*a + b*b )"

Changing a field using an expression

We create a 2D vector field with 2 values (vectors) and then we transform this field into a 3D vector field by applying an expression to values of the 2D field using applyFunc(). The expression func is applied to components of each vector of the field. We want the field to have 3 components computed as follows. (In func, we refer to the first component of a field value using the variable "a", and to the second component, using the variable "b").

• Component #0 = the second vector component hence "IVec * b" in func.
• Component #1 = the first vector component hence "JVec * a".
• Component #2 = a vector magnitude hence "KVec * sqrt( a*a + b*b )".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  // create a 2D vector field
  const double values[4] = {1.,1., 2.,1.};
  MCAuto<DataArrayDouble> array = DataArrayDouble::New();
  array->useExternalArrayWithRWAccess( values, 2, 2 ); // 2 tuples per 2 components
  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS );
  field->setArray( array );
  // transform the field to a 3D vector field
  const char func[] = "IVec * b + JVec * a + KVec * sqrt( a*a + b*b ) + 10";
  field->applyFunc( 3, func ); // require 3 components 
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 ); // 3 components as required

Now we ascertain that the result field is as we expect. We check the second vector of the field.

  double vec1[3]; // vector #1
  field->getArray()->getTuple( 1, vec1 );
  const double a = values[2], b = values[3]; // initial components of the vector #1
  CPPUNIT_ASSERT_DOUBLES_EQUAL( vec1[0], 10 + b, 13 ); // "10 + IVec * b"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( vec1[1], 10 + a, 13 ); // "10 + JVec * a"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( vec1[2], 10 + sqrt(a*a+b*b), 13 ); // "10 + KVec * sqrt( a*a + b*b )"

Filling a field with a value

We want to transform a 2D vector field to a 3D vector field so that all values to be equal to a certain value. First, we create the 2D mesh and the vector field on it.

Finally we use applyFunc() to change the number of components and all field values.

As a result, number of tuples in the field equals to the number of cells in the mesh, and number of components becomes equal to 3 as required.

Filling a field using an expression

First, we create a 2D Cartesian mesh constituted by 2 cells.

  const double coords[4] = {0.,2.,4.,6.}; // 6. is not used
  MCAuto<DataArrayDouble> x = DataArrayDouble::New();
  x->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<DataArrayDouble> y = DataArrayDouble::New();
  y->useExternalArrayWithRWAccess( coords, 2, 1 );
  MCAuto<MEDCouplingCMesh> mesh=MEDCouplingCMesh::New();
  mesh->setCoords(x,y);

Now we create a field on cells and use fillFromAnalytic2() to fill it with values computed using an expression func. This expression is applied to coordinates of each point (barycenter) for which the field value is computed. We want the field to have 3 components computed as follows. (In func, we refer to the first component of a point using the variable "a", and to the second component, using the variable "b").

• Component #0 = the second coordinate of the point hence "IVec * b" in func.
• Component #1 = the first coordinate of the point hence "JVec * a".
• Component #2 = distance between the point and SC origin (0.,0.) hence "KVec * sqrt( a*a + b*b )".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS );
  field->setMesh( mesh );
  const char func[] = "IVec * b + JVec * a + KVec * sqrt( a*a + b*b ) + 10";
  const char* varNames[2] = { "a", "b" }; // names used to refer to X and Y coord components
  std::vector<std::string> varNamesVec( varNames, varNames+2 );
  field->fillFromAnalyticNamedCompo( 3, varNamesVec, func );

Now we ascertain that the result field is as we expect. We check the second tuple of the field. We get barycenter of the cell #1 and checks that values of the second tuple are computed as we want.

  double val1[3]; // a value (vector) of the cell #1
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 ); // 3 components in the field
  field->getArray()->getTuple( 1, val1 );
  //
  MCAuto<DataArrayDouble> bc =
    mesh->computeCellCenterOfMass(); // func is applied to barycenters of cells
  double bc1[2]; // coordinates of the second point
  bc->getTuple( 1, bc1 );
  //
  double dist = sqrt( bc1[0]*bc1[0] + bc1[1]*bc1[1] );  // "sqrt( a*a + b*b )"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[0], 10 + bc1[1], 13 ); // "10 + IVec * b"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[1], 10 + bc1[0], 13 ); // "10 + JVec * a"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[2], 10 + dist  , 13 ); // "10 + KVec * sqrt( a*a + b*b )"

Filling a field using an expression

First, we create a 2D Cartesian mesh constituted by 2 cells. Note that we set names to coordinates arrays ("a" and "b" ) which will be used to refer to corresponding coordinates within a function.

  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> x = DataArrayDouble::New();
  x->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<DataArrayDouble> y = DataArrayDouble::New();
  y->useExternalArrayWithRWAccess( coords, 2, 1 );
  x->setInfoOnComponent(0,"a"); //  name used to refer to X coordinate within a function
  y->setInfoOnComponent(0,"b"); //  name used to refer to Y coordinate within a function
  MCAuto<MEDCouplingCMesh> mesh=MEDCouplingCMesh::New();
  mesh->setCoords(x,y);

Now we create a field on cells and use fillFromAnalytic2() to fill it with values computed using an expression func. This expression is applied to coordinates of each point (barycenter) for which the field value is computed. We want the field to have 3 components computed as follows. (In func, we refer to the first component of a point using the variable "a", and to the second component, using the variable "b").

• Component #0 = the second coordinate of the point hence "IVec * b" in func.
• Component #1 = the first coordinate of the point hence "JVec * a".
• Component #2 = distance between the point and SC origin (0.,0.) hence "KVec * sqrt( a*a + b*b )".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS );
  field->setMesh( mesh );
  const char func[] = "IVec * b + JVec * a + KVec * sqrt( a*a + b*b ) + 10";
  field->fillFromAnalytic( 3, func );

Now we ascertain that the result field is as we expect. We check the second tuple of the field. We get barycenter of the cell #1 and checks that values of the second tuple are computed as we want.

  double val1[3]; // a value (vector) of the cell #1
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 ); // 3 components in the field
  field->getArray()->getTuple( 1, val1 );
  //
  MCAuto<DataArrayDouble> bc =
    mesh->computeCellCenterOfMass(); // func is applied to barycenters of cells
  double bc1[2]; // coordinates of the second point
  bc->getTuple( 1, bc1 );
  //
  double dist = sqrt( bc1[0]*bc1[0] + bc1[1]*bc1[1] );  // "sqrt( a*a + b*b )"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[0], 10 + bc1[1], 13 ); // "10 + IVec * b"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[1], 10 + bc1[0], 13 ); // "10 + JVec * a"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[2], 10 + dist  , 13 ); // "10 + KVec * sqrt( a*a + b*b )"

Filling a field using an expression

First, we create a 2D Cartesian mesh constituted by 2 cells.

  const double coords[3] = {0.,2.,4};
  MCAuto<DataArrayDouble> x = DataArrayDouble::New();
  x->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<DataArrayDouble> y = DataArrayDouble::New();
  y->useExternalArrayWithRWAccess( coords, 2, 1 );
  MCAuto<MEDCouplingCMesh> mesh=MEDCouplingCMesh::New();
  mesh->setCoords(x,y);

Now we create a field on cells and use fillFromAnalytic() to fill it with values computed using an expression func. This expression is applied to coordinates of each point (barycenter) for which the field value is computed. We want the field to have 3 components computed as follows. (In func, we refer to the first component of a point using the variable "a", and to the second component, using the variable "b").

• Component #0 = the second coordinate of the point hence "IVec * b" in func.
• Component #1 = the first coordinate of the point hence "JVec * a".
• Component #2 = distance between the point and SC origin (0.,0.) hence "KVec * sqrt( a*a + b*b )".

In addition we want to add 10.0 to each component computed as described above, hence "10" in func.

  const char func[] = "IVec * b + JVec * a + KVec * sqrt( a*a + b*b ) + 10";
  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS );
  field->setMesh( mesh );
  field->fillFromAnalytic( 3, func );

Now we ascertain that the result field is as we expect. We check the second tuple of the field. We get barycenter of the cell #1 to check that values of the second tuple (#1) are computed as we want.

  double val1[3]; // a value (vector) of the cell #1
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 ); // 3 components in the field
  field->getArray()->getTuple( 1, val1 );
  //
  MCAuto<DataArrayDouble> bc =
    mesh->computeCellCenterOfMass(); // func is applied to barycenters of cells
  double bc1[2]; // coordinates of the second point
  bc->getTuple( 1, bc1 );
  //
  double dist = sqrt( bc1[0]*bc1[0] + bc1[1]*bc1[1] );  // "sqrt( a*a + b*b )"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[0], 10 + bc1[1], 13 ); // "10 + IVec * b"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[1], 10 + bc1[0], 13 ); // "10 + JVec * a"
  CPPUNIT_ASSERT_DOUBLES_EQUAL( val1[2], 10 + dist  , 13 ); // "10 + KVec * sqrt( a*a + b*b )"

Some operations that can be carried out on fields on cells

  MEDCoupling::MEDCouplingFieldDouble *f1=mesh->fillFromAnalytic(MEDCoupling::ON_CELLS,1,"x*x+y*y*3+2.*x");//f1 is scalar
  MEDCoupling::MEDCouplingFieldDouble *f2=mesh->fillFromAnalytic(MEDCoupling::ON_CELLS,1,"cos(x+y/x)");//f2 is scalar too
  MEDCoupling::MEDCouplingFieldDouble *f2bis=mesh->fillFromAnalytic(MEDCoupling::ON_CELLS,2,"x*x*IVec+3*y*JVec");//f2bis is a vectors field
  MEDCoupling::MEDCouplingFieldDouble *f3=(*f1)+(*f2);//f3 scalar
  MEDCoupling::MEDCouplingFieldDouble *f4=(*f3)/(*f2);//f4 scalar
  f2bis->applyFunc(1,"sqrt(x*x+y*y)");//f2bis becomes scalar
  MEDCoupling::MEDCouplingFieldDouble *f5=(*f2bis)*(*f4);//f5 scalar
  const double pos1[2]={0.48,0.38};
  double res;
  f4->getValueOn(pos1,&res);//f4 is scalar so the returned value is of size 1.
  // ...

The decrementation of ref counter should be carried out in CPlusPlus only ...

  // f1, f2, f2bis, f3, f4, f5 are no more useful here : release them
  f1->decrRef();
  f2->decrRef();
  f2bis->decrRef();
  f3->decrRef();
  f4->decrRef();
  f5->decrRef();

Permuting a field on nodes

First, we create a supporting 2D mesh constituted by 4 cells. We create a 2x2 Cartesian mesh and then convert it to an unstructured one, since the Cartesian mesh is not suitable for renumberNodes() as its nature does not imply node renumbering.

  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> cmesh = MEDCouplingCMesh::New();
  cmesh->setCoords(coordsArr,coordsArr);
  MCAuto<MEDCouplingUMesh> mesh = cmesh->buildUnstructured();

Then we create a field on nodes using fillFromAnalytic(), such that its values to coincide with coordinates of field location points that are nodes in our case (as our field is ON_NODES). At last we ascertain that field values are equal to node coordinates.

  MCAuto<MEDCouplingFieldDouble> field =
    mesh->fillFromAnalytic( MEDCoupling::ON_NODES,2,"IVec*x+JVec*y");
  const DataArrayDouble* values = field->getArray();
  const DataArrayDouble* nodeCoords = mesh->getCoords();
  CPPUNIT_ASSERT( values->isEqualWithoutConsideringStr( *nodeCoords, 1e-13 ));

Now, we are going to reverse order of nodes using renumberNodes().

  const mcIdType renumber[9] = { 8, 7, 6, 5, 4, 3, 2, 1, 0 };
  field->renumberNodes(renumber,false);
  const MEDCouplingMesh* mesh2 = field->getMesh(); // field now refers to another mesh
  values = field->getArray();
  nodeCoords = (static_cast<const MEDCouplingUMesh*>(mesh2))->getCoords();
  CPPUNIT_ASSERT( values->isEqualWithoutConsideringStr( *nodeCoords, 1e-13 ));

As a result, the underlying mesh of field is changed and its nodes are also renumbered. And the field values are still equal to node coordinates of the renumbered mesh2.

Permuting a field on cells

First, we create a supporting 2D mesh constituted by 4 cells. We create a 2x2 Cartesian mesh and then convert it to an unstructured one, since the Cartesian mesh is not suitable for renumberCells() as its nature does not imply cell renumbering.

  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> cmesh = MEDCouplingCMesh::New();
  cmesh->setCoords(coordsArr,coordsArr);
  MCAuto<MEDCouplingUMesh> mesh = cmesh->buildUnstructured();

Then we create a field on cells using fillFromAnalytic(), such that its values to coincide with coordinates of field location points that are cell barycenters in our case (as our field is ON_CELLS). At last we ascertain that field values are equal to cell barycenters.

  MCAuto<MEDCouplingFieldDouble> field =
    mesh->fillFromAnalytic( MEDCoupling::ON_CELLS,2,"IVec*x+JVec*y");
  const DataArrayDouble* values = field->getArray();
  MCAuto<DataArrayDouble> bc = mesh->computeCellCenterOfMass();
  CPPUNIT_ASSERT( values->isEqualWithoutConsideringStr( *bc, 1e-13 ));

Now, we are going to reverse order of cells using renumberCells().

  const mcIdType renumber[4] = { 3, 2, 1, 0 };
  field->renumberCells(renumber,false);
  const MEDCouplingMesh* mesh2 = field->getMesh(); // field now refers to another mesh
  values = field->getArray();
  bc = mesh2->computeCellCenterOfMass();
  CPPUNIT_ASSERT( values->isEqualWithoutConsideringStr( *bc, 1e-13 ));

As a result, the underlying mesh of field is changed and its cells are also renumbered. And the field values are still equal to cell barycenters of the renumbered mesh2.

Filling a field using a C function

We want to create a 3D vector field lying on a 2D mesh using a C function as a value generator. For that, first, we define the function that computes 3 values basing on 2 coordinates of a 2D point.

bool getNewValue(const double *pos, double *res)
{
  res[0] = pos[0];
  res[1] = pos[1];
  res[2] = sqrt( pos[0]*pos[0] + pos[1]*pos[1] );
  return true;
}

Then we create the 2D mesh and the field on it, and finally we use fillFromAnalytic() to fill the field with values each composed of 3 components.

  // mesh
  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh = MEDCouplingCMesh::New();
  mesh->setCoords(coordsArr,coordsArr); // mesh becomes a 2D structured mesh
  // field
  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS );
  field->setMesh( mesh );
  field->fillFromAnalytic( 3, &getNewValue ); // 3 components are required
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 );
  CPPUNIT_ASSERT( field->getNumberOfTuples() == mesh->getNumberOfCells() );

As a result, number of tuples in the field equals to the number of cells in the mesh, and number of components equals to 3 as required.

Changing a field by applying a C function

We want to transform a 2D vector field to a 3D vector field by applying a C function to each vector value. For that, first, we define the function that computes 3 values basing on 2 components of a 2D vector.

bool getNewValue(const double *pos, double *res)
{
  res[0] = pos[0];
  res[1] = pos[1];
  res[2] = sqrt( pos[0]*pos[0] + pos[1]*pos[1] );
  return true;
}

Then we create the 2D mesh and the vector field on it.

  // mesh
  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh = MEDCouplingCMesh::New();
  mesh->setCoords(coordsArr,coordsArr); // mesh becomes a 2D structured mesh
  // field
  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS );
  field->setMesh( mesh );
  MCAuto<DataArrayDouble> bc = mesh->computeCellCenterOfMass();
  field->setArray( bc ); // 2 components here as the mesh is 2D

Finally we use applyFunc() to change the field values.

  field->applyFunc( 3, &getNewValue ); // 3 components are required
  CPPUNIT_ASSERT( field->getNumberOfComponents() == 3 );
  CPPUNIT_ASSERT( field->getNumberOfTuples() == mesh->getNumberOfCells() );

As a result, number of tuples in the field equals to the number of cells in the mesh, and number of components becomes equal to 3 as required.

Using interpolation tools - High level usage

...
const char sourceFileName[]="source.med"
MEDCouplingFieldDouble *sourceField=MEDLoader::ReadFieldCell(sourceFileName,"Source_Mesh",0,"Density",/*iteration*/0,/*order*/0)
const char targetFileName[]="target.med"
MEDCouplingUMesh *med_target_mesh=MEDLoader::ReadUMeshFromFile(targetFileName,"Target_Mesh",0)
//
sourceField->setNature(ConservativeVolumic)//Specify which formula to use in case of non overlapping meshes
MEDCouplingRemapper remapper
remapper.setPrecision(1e-12)
remapper.setIntersectionType(INTERP_KERNEL::Triangulation)
remapper.prepare(sourceField->getMesh(),med_target_mesh,"P0P0")
MEDCouplingFieldDouble *targetField=remapper.transferField(sourceField,/*default_value*/4.57)//Any target cell not intercepted by any source cell will have value set to 4.57.
...
// clean-up
targetField->decrRef()
sourceField->decrRef()
med_target_mesh->decrRef()

Using interpolation tools - Middle level usage

• The simplest way to use the interpolator with MEDCoupling data structure is illustrated in the following example.

...
MEDCouplingUMesh *med_source_mesh=MEDLoader::ReadUMeshFromFile("source.med","Source_mesh",0)
MEDCouplingUMesh *med_target_mesh=MEDLoader::ReadUMeshFromFile("target.med","Target_mesh",0)
MEDCouplingNormalizedUnstructuredMesh<2,2> wrap_source_mesh(med_source_mesh)
MEDCouplingNormalizedUnstructuredMesh<2,2> wrap_target_mesh(med_target_mesh)
// Go for interpolation...
INTERP_KERNEL::Interpolation2D myInterpolator
myInterpolator.setPrecision(1e-7)
myInterpolator.setIntersectionType(INTERP_KERNEL::Geometric2D)
std::vector<std::map<int,double> > resultMatrix
INTERP_KERNEL::Matrix<double,ALL_C_MODE> resultMatrix2
// here the interpolation is performed twice for this code to illustrate the possibility of storing data the interpolation matrix in 2 different data structures.
myInterpolator.interpolateMeshes(wrap_source_mesh,wrap_target_mesh,resultMatrix,"P0P0")
myInterpolator.interpolateMeshes(wrap_source_mesh,wrap_target_mesh,resultMatrix2,"P0P0")
//Ok resultMatrix and resultMatrix2 contain matrix now
...

• Same with VTK data structure :

...
vtkXMLUnstructuredGridReader *readerSource=vtkXMLUnstructuredGridReader::New()
readerSource->SetFileName("source.vtu")
vtkUnstructuredGrid *vtk_source_mesh=readerSource->GetOutput()
readerSource->Update()
vtkXMLUnstructuredGridReader *readerTarget=vtkXMLUnstructuredGridReader::New()
readerTarget->SetFileName("target.vtu")
vtkUnstructuredGrid *vtk_target_mesh=readerTarget->GetOutput()
readerTarget->Update()
// Ok at this point we have our mesh in VTK format.
// Go to wrap vtk_source_mesh and vtk_target_mesh.
VTKNormalizedUnstructuredMesh<2> wrap_source_mesh(vtk_source_mesh)
VTKNormalizedUnstructuredMesh<2> wrap_target_mesh(vtk_target_mesh)
// Go for interpolation...
INTERP_KERNEL::Interpolation2D myInterpolator
//optional call to parametrize your interpolation. First precision, tracelevel, intersector wanted.
myInterpolator.setOptions(1e-7,0,Geometric2D)
INTERP_KERNEL::Matrix<double,ALL_C_MODE> resultMatrix
myInterpolator.interpolateMeshes(wrap_source_mesh,wrap_target_mesh,resultMatrix,"P0P0")
//Ok let's multiply resultMatrix by source field to interpolate to target field.
resultMatrix.multiply(...)
//clean-up
readerSource->Delete()
readerTarget->Delete()
...

Access

Getting a field value at some point at certain time

First, we create a supporting structured mesh. We create a 2x2 Cartesian mesh constituted by 4 cells.

  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh = MEDCouplingCMesh::New();
  mesh->setCoords(coordsArr,coordsArr);

Then we create a scalar field on cells, whose values vary linearly in time. We set all field values at a start time to be equal 10.0 using fillFromAnalytic(). And we set all field values at an end time to be equal 20.0 by doubling the start time array.

  MCAuto<MEDCouplingFieldDouble> field =
    MEDCouplingFieldDouble::New( MEDCoupling::ON_CELLS, MEDCoupling::LINEAR_TIME );
  field->setMesh( mesh );
  field->fillFromAnalytic( 1,"10"); // all values == 10.
  MCAuto<DataArrayDouble> array2 =
    DataArrayDouble::Add( field->getArray(), field->getArray() ); // == 2 * field->getArray()
  field->setEndArray( array2 ); // all values == 20.
  const double time1 = 1.1, time2 = 22.;
  field->setStartTime( time1, 0, 0 );
  field->setEndTime  ( time2, 0, 0 );

Now, we want to get a field value at a point [0,0] at a middle time between the start and end times. We expect the returned value to be equal to an average of 10. and 20.

  const double pos[2] = { 1., 1. }; // we are in 2D space
  double value[1]; // the field is scalar <-> 1 component
  field->getValueOn( pos, 0.5*( time1 + time2 ), value );
  CPPUNIT_ASSERT( fabs( value[0] - 0.5*( 10. + 20. )) < 1e-13 ); 

Getting field values at some points

First, we create a supporting structured mesh. We create a 2x2 Cartesian mesh constituted by 4 cells. Then we create a scalar field on cells using fillFromAnalytic().

  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh = MEDCouplingCMesh::New();
  mesh->setCoords(coordsArr,coordsArr);
  MCAuto<MEDCouplingFieldDouble> field =
    mesh->fillFromAnalytic( MEDCoupling::ON_CELLS,1,"x+y");

Now, we want to retrieve all field values using getValueOnMulti(). The field values relate to cells, hence we will use cell barycenters as a parameter of getValueOnMulti(). We expect that the double array returned getValueOnMulti() is equal to that stored by field.

  // field values are located at cell barycenters
  MCAuto<DataArrayDouble> bc = mesh->computeCellCenterOfMass();
  MCAuto<DataArrayDouble> valArray =
    field->getValueOnMulti( bc->getConstPointer(), bc->getNumberOfTuples() );
  CPPUNIT_ASSERT( valArray->isEqual( * field->getArray(), 1e-13 ));

Getting a field value at a point

First, we create a supporting structured mesh. We create a 2x2 Cartesian mesh constituted by 4 cells. Then we create a scalar field on cells using fillFromAnalytic().

  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh = MEDCouplingCMesh::New();
  mesh->setCoords(coordsArr,coordsArr);
  MCAuto<MEDCouplingFieldDouble> field =
    mesh->fillFromAnalytic( MEDCoupling::ON_CELLS,1,"x+y");

Now, we want to retrieve all field values using getValueOn(). The field values relate to cells, hence we will use cell barycenters to get a field value at each cell.

  // field values are located at cell barycenters
  MCAuto<DataArrayDouble> bc = mesh->computeCellCenterOfMass();
  std::vector<double> vals( field->getNumberOfTuples() ); // array to collect values returned by getValueOn()
  double cellBC[2]; // we are in 2D space
  for ( int i = 0; i < bc->getNumberOfTuples(); ++i )
  {
    bc->getTuple( i, cellBC );
    field->getValueOn( cellBC, & vals[i] );
  }
  CPPUNIT_ASSERT( std::equal( vals.begin(), vals.end(), field->getArray()->getConstPointer() ));

We collected all values returned by getValueOn() in an array, so that we can ascertain that the array of returned values is same as that stored by field.

Getting a value of field lying on a structured mesh

First, we create a supporting structured mesh. We create a 2x2 Cartesian mesh constituted by 4 cells. Then we create a scalar field on cells using fillFromAnalytic().

  const double coords[4] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh = MEDCouplingCMesh::New();
  mesh->setCoords(coordsArr,coordsArr);
  MCAuto<MEDCouplingFieldDouble> field =
    mesh->fillFromAnalytic( MEDCoupling::ON_CELLS,1,"x+y");

Now, we retrieve a field value relating to the cell #3 (this cell has a structured indexed (1,1)). For that we use getValueOnPos() where we pass the structured indexed of the cell: 1,1,-1 (the last index is meaningless as the mesh is 2D).

  double val11[1]; // 1 == field->getNumberOfComponents()
  field->getValueOnPos( 1,1,-1, val11 );
  // field values are located at cell barycenters
  MCAuto<DataArrayDouble> bc = mesh->computeCellCenterOfMass();
  CPPUNIT_ASSERT( val11[0] == bc->getIJ(3,0) + bc->getIJ(3,1) );

After all we ascertain that the returned value corresponds to the expression used for the field creation. Namely that the value equals to the sum of components of barycenter of cell #3.

Meshes

Create

Standard build of an unstructured mesh from scratch

Firstly retrieve basic data in full interlace mode for coordinates, and nodal connectivity cell per cell.

  double coords[27]={-0.3,-0.3,0.,   0.2,-0.3,0.,   0.7,-0.3,0.,   -0.3,0.2,0.,   0.2,0.2,0., 
                     0.7,0.2,0.,    -0.3,0.7,0.,    0.2,0.7,0.,     0.7,0.7,0. };
  mcIdType nodalConnPerCell[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};

Then create MEDCoupling::MEDCouplingUMesh instance giving its mesh dimension (2 here) and a name.

  MEDCoupling::MEDCouplingUMesh *mesh=MEDCoupling::MEDCouplingUMesh::New("My2DMesh",2);

Gives an upper bound of the number of cells to be inserted into the unstructured mesh.
Then enter nodal connectivity of all cells, cell per cell using MEDCoupling::MEDCouplingUMesh::insertNextCell method.
When the nodal connectivity cell per cell has been finished, call MEDCoupling::MEDCouplingUMesh::finishInsertingCells method in order to restore mesh instance.

  mesh->allocateCells(5);//You can put more than 5 if you want but not less.
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,nodalConnPerCell);
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,nodalConnPerCell+4);
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,nodalConnPerCell+7);
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,nodalConnPerCell+10);
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,nodalConnPerCell+14);
  mesh->finishInsertingCells();

At this level the connectivity part of the mesh mesh as been defined. Now let's set the coordinates using array coords defined above.

  MEDCoupling::DataArrayDouble *coordsArr=MEDCoupling::DataArrayDouble::New();
  coordsArr->alloc(9,3);//here coordsArr are declared to have 3 components, mesh will deduce that its spaceDim==3. 
  std::copy(coords,coords+27,coordsArr->getPointer());
  mesh->setCoords(coordsArr);//coordsArr contains 9 tuples, that is to say mesh contains 9 nodes.
  coordsArr->decrRef();

At this level mesh is usable. When this mesh is no more needed simply call decrRef to decrement its reference counter.

  mesh->decrRef();

Advanced build of an unstructured mesh from scratch

Firstly retrieve basic data in full interlace mode for coordinates, and nodal connectivity cell per cell, cell type included (3 for INTERP_KERNEL::NORM_TRI3 and 4 for INTERP_KERNEL::QUAD4).

  double coords[27]={-0.3,-0.3,0.,   0.2,-0.3,0.,   0.7,-0.3,0.,   -0.3,0.2,0.,   0.2,0.2,0., 
                     0.7,0.2,0.,    -0.3,0.7,0.,    0.2,0.7,0.,     0.7,0.7,0. };
  mcIdType nodalConnPerCell[23]={4,0,3,4,1, 3,1,4,2, 3,4,5,2, 4,6,7,4,3, 4,7,8,5,4};
  mcIdType nodalConnPerCellIndex[6]={0,5,9,13,18,23};

Then create MEDCoupling::MEDCouplingUMesh instance giving its mesh dimension (2 here) and a name.

  MEDCoupling::MEDCouplingUMesh *mesh=MEDCoupling::MEDCouplingUMesh::New("My2DMesh",2);

Then enter nodal connectivity at once.

  MEDCoupling::DataArrayIdType *nodalConn=MEDCoupling::DataArrayIdType::New();
  nodalConn->alloc(23,1);
  std::copy(nodalConnPerCell,nodalConnPerCell+23,nodalConn->getPointer());
  MEDCoupling::DataArrayIdType *nodalConnI=MEDCoupling::DataArrayIdType::New();
  nodalConnI->alloc(6,1);
  std::copy(nodalConnPerCellIndex,nodalConnPerCellIndex+6,nodalConnI->getPointer());
  mesh->setConnectivity(nodalConn,nodalConnI,true);
  nodalConn->decrRef();// nodalConn DataArrayIdType instance is owned by mesh after call to setConnectivity method. No more need here -> decrRef()
  nodalConnI->decrRef();// nodalConnI DataArrayIdType instance is owned by mesh after call to setConnectivity method. No more need here -> decrRef()

At this level the connectivity part of the mesh mesh as been defined. Now let's set the coordinates using array coords defined above.

  MEDCoupling::DataArrayDouble *coordsArr=MEDCoupling::DataArrayDouble::New();
  coordsArr->alloc(9,3);//here coordsArr are declared to have 3 components, mesh will deduce that its spaceDim==3. 
  std::copy(coords,coords+27,coordsArr->getPointer());
  mesh->setCoords(coordsArr);//coordsArr contains 9 tuples, that is to say mesh contains 9 nodes.
  coordsArr->decrRef();

At this level mesh is usable. When this mesh is no more needed simply call decrRef() to decrement its reference counter.

  mesh->decrRef();

Standard build of an cartesian mesh from scratch

We are going to build a 2D cartesian mesh, constituted from 9 nodes along X axis, and 7 nodes along Y axis.

Firstly retrieve for each direction the discretization and build a DataArrayDouble instance on the corresponding direction.

  double XCoords[9]={-0.3,0.,0.1,0.3,0.45,0.47,0.49,1.,1.22};
  double YCoords[7]={0.,0.1,0.37,0.45,0.47,0.49,1.007};
  MEDCoupling::DataArrayDouble *arrX=MEDCoupling::DataArrayDouble::New();
  arrX->alloc(9,1);
  std::copy(XCoords,XCoords+9,arrX->getPointer());
  arrX->setInfoOnComponent(0,"X [m]");
  MEDCoupling::DataArrayDouble *arrY=MEDCoupling::DataArrayDouble::New();
  arrY->alloc(7,1);
  std::copy(YCoords,YCoords+7,arrY->getPointer());
  arrY->setInfoOnComponent(0,"Y [m]");

Then create MEDCoupling::MEDCouplingCMesh instance giving the 2 instances of DataArrayDouble obtained above.

There are 2 techniques to get it.

Either :

  MEDCoupling::MEDCouplingCMesh *mesh=MEDCoupling::MEDCouplingCMesh::New("My2D_CMesh");
  mesh->setCoords(arrX,arrY);
  arrX->decrRef();
  arrY->decrRef();

Or :

  mesh->setCoordsAt(0,arrX);
  arrX->decrRef();
  mesh->setCoordsAt(1,arrY);
  arrY->decrRef();

mesh is now available for use :

  CPPUNIT_ASSERT_EQUAL(ToIdType(8*6),mesh->getNumberOfCells());
  CPPUNIT_ASSERT_EQUAL(ToIdType(9*7),mesh->getNumberOfNodes());
  CPPUNIT_ASSERT_EQUAL(2,mesh->getSpaceDimension());
  CPPUNIT_ASSERT_EQUAL(2,mesh->getMeshDimension());

When this mesh is no more needed simply call decrRef to decrement its reference counter (nothing to be done in Python).

  mesh->decrRef();

Getting a bounding mesh

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);   
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10);
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14);
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we use buildBoundaryMesh() to get a mesh of lower dimension bounding mesh.

  MCAuto<MEDCouplingPointSet> mesh1=mesh->buildBoundaryMesh(true);
  MCAuto<MEDCouplingPointSet> mesh2=mesh->buildBoundaryMesh(false);
  CPPUNIT_ASSERT(  coordsArr->isEqual( *mesh1->getCoords(), 1e-13 )); // same nodes
  CPPUNIT_ASSERT( !coordsArr->isEqual( *mesh2->getCoords(), 1e-13 )); // different nodes

Depending on the value of a parameter, buildBoundaryMesh() creates the mesh sharing the node coordinates array with mesh or not.

Retrieving a lower dimension mesh based on given nodes

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

In the following code we retrieve nodes of the cell #0 an then we call buildFacePartOfMySelfNode() twice with these nodes and with varying last parameter allNodes as input.

  std::vector<mcIdType> nodes;
  mesh->getNodeIdsOfCell( 0, nodes );
  const bool allNodes = true;
  MCAuto<MEDCouplingUMesh> mesh1 =
    (MEDCouplingUMesh*)mesh->buildFacePartOfMySelfNode( &nodes[0],&nodes[0]+nodes.size(),allNodes);
  CPPUNIT_ASSERT( mesh1->getNumberOfCells() == 4 ); // 4 segments bounding QUAD4 #0 only
  MCAuto<MEDCouplingUMesh> mesh2 =
    (MEDCouplingUMesh*)mesh->buildFacePartOfMySelfNode( &nodes[0],&nodes[0]+nodes.size(),!allNodes);
  CPPUNIT_ASSERT( mesh2->getNumberOfCells() == 9 ); // more segments added

If the last parameter is true buildFacePartOfMySelfNode() looks for segments whose all nodes are given to it, hence it finds segments bounding the cell #0 only.
If the last parameter is false buildFacePartOfMySelfNode() looks for any segment whose nodes are given to it, hence it adds more segments to mesh2.

Copying cells selected by nodes

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

In the following code we retrieve nodes of the cell #0 an then we call buildPartOfMySelfNode() twice with these nodes and with varying last parameter allNodes as input.

  std::vector<mcIdType> nodes;
  mesh->getNodeIdsOfCell( 0, nodes );
  const bool allNodes = true;
  MCAuto<MEDCouplingUMesh> mesh1 =
    (MEDCouplingUMesh*)mesh->buildPartOfMySelfNode( &nodes[0], &nodes[0]+nodes.size(), allNodes);
  MCAuto<MEDCouplingUMesh> mesh2 =
    (MEDCouplingUMesh*)mesh->buildPartOfMySelfNode( &nodes[0], &nodes[0]+nodes.size(),!allNodes);
  CPPUNIT_ASSERT_EQUAL( mesh1->getNumberOfCells(), ToIdType( 1 ));
  CPPUNIT_ASSERT_EQUAL( mesh2->getNumberOfCells(), mesh->getNumberOfCells() );

If the last parameter is true buildPartOfMySelfNode() looks for cells whose all nodes are given to it, hence it finds the cell #0 only.
If the last parameter is false buildPartOfMySelfNode() looks for any cell whose nodes are given to it, hence it finds all cells of mesh because all cells share the node #4.

Getting a part of mesh

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we use buildPartOfMySelf() to get a mesh containing only two cells of mesh.

  const mcIdType cellIds[2]={1,2};
  MEDCouplingUMesh* mesh2=(MEDCouplingUMesh*)mesh->buildPartOfMySelf(cellIds,cellIds+2,true);
  MEDCouplingUMesh* mesh3=(MEDCouplingUMesh*)mesh->buildPartOfMySelf(cellIds,cellIds+2,false);
  CPPUNIT_ASSERT(  coordsArr->isEqual( *mesh2->getCoords(), 1e-13 )); // same nodes
  CPPUNIT_ASSERT( !coordsArr->isEqual( *mesh3->getCoords(), 1e-13 )); // different nodes
  for ( mcIdType i = 0; i < 2; ++i )
    {
      std::vector<mcIdType> nodes1, nodes2;
      mesh ->getNodeIdsOfCell(cellIds[i], nodes1);
      mesh2->getNodeIdsOfCell(i, nodes2);
      CPPUNIT_ASSERT( nodes1 == nodes2 ); // cell #cellIds[i] was copied
    }

Modify

Fixing orientation of "extruded" volumes

First, we create a mesh with 2 incorrectly oriented "extruded" volumes.

  // 2D coordinates of 5 base nodes
  const double coords[5*2]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2 };
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 5, 2 );
  // coordinates of 5 top nodes
  MCAuto<DataArrayDouble> coordsArr2 = coordsArr->deepCopy();
  // 3D coordinates of base + top nodes
  coordsArr  = coordsArr-> changeNbOfComponents( 3, 0 );
  coordsArr2 = coordsArr2->changeNbOfComponents( 3, 1 );
  coordsArr = DataArrayDouble::Aggregate( coordsArr, coordsArr2 );
  // mesh
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);
  mesh->setMeshDimension(3);
  mesh->allocateCells(2);
  // connectivity of reversed HEXA8 and PENTA6
  const mcIdType conn[8+6]={0,1,4,3, 5,6,9,8, 1,2,4, 6,7,9};
  mesh->insertNextCell(INTERP_KERNEL::NORM_HEXA8, 8,conn+0);
  mesh->insertNextCell(INTERP_KERNEL::NORM_PENTA6,6,conn+8);
  mesh->finishInsertingCells();

Now we check that findAndCorrectBadOriented3DExtrudedCells() finds and fixes the reversed cells.

  MCAuto<DataArrayIdType> fixedCells =
    mesh->findAndCorrectBadOriented3DExtrudedCells();
  CPPUNIT_ASSERT( fixedCells->getNumberOfTuples() == 2 ); // 2 cells fixed
  fixedCells = mesh->findAndCorrectBadOriented3DExtrudedCells();
  CPPUNIT_ASSERT( fixedCells->getNumberOfTuples() == 0 ); // no bad cells

Fixing orientation of polyhedra

First, we create a mesh with 2 polyhedra, one of which is incorrectly oriented. We create two "extruded" polyhedra and then convert them to correctly defined polyhedra.

  // 2D coordinates of 5 base nodes
  const double coords[5*2]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2 };
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 5, 2 );
  // coordinates of 5 top nodes
  MCAuto<DataArrayDouble> coordsArr2 = coordsArr->deepCopy();
  // 3D coordinates of base + top nodes
  coordsArr  = coordsArr-> changeNbOfComponents( 3, 0 );
  coordsArr2 = coordsArr2->changeNbOfComponents( 3, 1 );
  coordsArr = DataArrayDouble::Aggregate( coordsArr, coordsArr2 );
  // mesh
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);
  mesh->setMeshDimension(3);
  mesh->allocateCells(2);
  // connectivity of a HEXA8 + a reversed PENTA6
  const mcIdType conn[8+6]={0,3,4,1, 5,8,9,6, 1,2,4, 6,7,9};
  mesh->insertNextCell(INTERP_KERNEL::NORM_POLYHED,8,conn); //  "extruded" polyhedron
  mesh->insertNextCell(INTERP_KERNEL::NORM_POLYHED,6,conn+8);
  mesh->finishInsertingCells();
  // fix connectivity of NORM_POLYHED's
  mesh->convertExtrudedPolyhedra();

Now we check that arePolyhedronsNotCorrectlyOriented() finds one reversed cell. After that we fix it using orientCorrectlyPolyhedrons() and re-check the orientation of polyhedra.

  std::vector<mcIdType> badCellIds;
  mesh->arePolyhedronsNotCorrectlyOriented( badCellIds );
  CPPUNIT_ASSERT( badCellIds.size() == 1 ); //  one polyhedron is KO
  // fix invalid rolyherdons
  mesh->orientCorrectlyPolyhedrons();
  // re-check orientation
  badCellIds.clear(); // as badCellIds is not cleared by arePolyhedronsNotCorrectlyOriented()
  mesh->arePolyhedronsNotCorrectlyOriented( badCellIds );
  CPPUNIT_ASSERT( badCellIds.size() == 0 ); // connectivity is OK

Fixing orientation of faces

First, we create a 2D mesh in 3D space with 3 QUAD4 and 2 TRI3 cells. Orientation of the cell #1 is reversed comparing with others.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,2,4, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);
  mesh->changeSpaceDimension(3);

Now we check that are2DCellsNotCorrectlyOriented() finds one reversed face. After that we fix the incorrectly oriented cell using orientCorrectly2DCells() and re-check the orientation of cells.

  const double vec[3] = {0.,0.,-1.};
  std::vector<mcIdType> badCellIds;
  mesh->are2DCellsNotCorrectlyOriented( vec, false, badCellIds );
  CPPUNIT_ASSERT( badCellIds.size() == 1 ); //  one cell is reversed
  // fix orientation
  mesh->orientCorrectly2DCells( vec, false );
  // re-check orientation
  badCellIds.clear(); // as badCellIds is not cleared by are2DCellsNotCorrectlyOriented()
  mesh->are2DCellsNotCorrectlyOriented( vec, false, badCellIds );
  CPPUNIT_ASSERT( badCellIds.size() == 0 ); // the orientation is OK

Alternatively you can orient all 2D cells equally using the first cell as a reference:

  mesh->orientCorrectly2DCells();

Also it is possible to orient some selected 2D cells by using another group of cells as the reference:

  const mcIdType refCells[] = { 0,2 };
  const mcIdType objCells[] = { 1,3 };
  MCAuto<MEDCouplingUMesh> refGroup = mesh->buildPartOfMySelf( refCells, refCells + 2 );
  MCAuto<MEDCouplingUMesh> objGroup = mesh->buildPartOfMySelf( objCells, objCells + 2 );
  objGroup->orientCorrectly2DCells( refGroup );
  mesh->setPartOfMySelf( objCells, objCells + 2, *objGroup );

Renumbering nodes in the connectivity array

First, we create a 2D mesh with 1 QUAD4 cell and with undefined coordinates of nodes.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(1);
  const mcIdType conn[4]={4,3,2,1};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);
  mesh->finishInsertingCells();

Now we use renumberNodesInConn() to get the following nodal connectivity of a sole cell: 0,1,2,3.

  const mcIdType old2newIds[] = {-1,3,2,1,0};
  mesh->renumberNodesInConn( old2newIds );
  const mcIdType nodes0Expected[] = {0,1,2,3};
  std::vector<mcIdType> nodes0;
  mesh->getNodeIdsOfCell( 0, nodes0 );
  CPPUNIT_ASSERT(std::equal( nodes0Expected, nodes0Expected+4, &nodes0[0] ));

old2newIds array defines how node ids are changed:

• new id of node #0 is -1,
• new id of node #1 is 3,
• new id of node #2 is 4,
• new id of node #3 is 1,
• new id of node #4 is 0.

Renumbering nodes

First, we create a 2D mesh with 4 nodes and no cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  const double coords[4*2]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.3};
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(coords, 4,2);
  mesh->setCoords(coordsArr);
  mesh->allocateCells(0);
  mesh->finishInsertingCells();

Next, we use renumberNodes() to permute nodes so that

• old node #0 becomes #2,
• old node #1 remains #1,
• old node #2 becomes #0,
• old node #3 is removed.

Number of nodes becomes 3.

  const mcIdType newIds[] = { 2,1,0,-1 };
  mesh->renumberNodes(newIds, 3);
  coordsArr = mesh->getCoordinatesAndOwner(); // get a shorten array
  const double coordsExpected[3*2]={0.7,-0.3, 0.2,-0.3, -0.3,-0.3};
  MCAuto<DataArrayDouble> coordsExpectedArr=DataArrayDouble::New();
  coordsExpectedArr->useExternalArrayWithRWAccess(coordsExpected, 3,2);
  CPPUNIT_ASSERT( coordsExpectedArr->isEqual( *coordsArr, 1e-13 ));

Next we compare behavior of renumberNodes() and that of renumberNodes2() which, in contrast to renumberNodes(), moves merged nodes to their barycenter.
We set #2 as new id of old node #3 and expect that renumberNodes2() moves old nodes #0 and #3 to their barycenter (-0.3,0.0) which becomes position of node #2.

  coordsArr->useExternalArrayWithRWAccess(coords, 4,2); // restore old nodes
  const mcIdType newIds2[] = { 2,1,0,2 };
  mesh->renumberNodesCenter(newIds2, 3);
  coordsArr = mesh->getCoordinatesAndOwner(); // get a shorten array
  const double coordsExpected2[3*2]={0.7,-0.3, 0.2,-0.3, -0.3, 0.0};
  MCAuto<DataArrayDouble> coordsExpectedArr2=DataArrayDouble::New();
  coordsExpectedArr2->useExternalArrayWithRWAccess(coordsExpected2, 3,2);
  CPPUNIT_ASSERT( coordsExpectedArr2->isEqual( *coordsArr, 1e-13 ));

Merging equal nodes

First, we create a 2D mesh with 1 QUAD4 and 2 TRI3 cells. The cells are based on 6 nodes of which 2 nodes fully coincide (#3 and #4) and 3 nodes are equal with precision 0.003.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);
  mesh->finishInsertingCells();
  const double coords[6*2]={0.3,-0.301,  // #0
                            0.2,-0.3,    // #1
                            0.3,-0.302,  // #2 ~~ #0
                            1.1,0.0,     // #3
                            1.1,0.0,     // #4 == #3
                            0.3,-0.303}; // #5 ~~ #0
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(6,2);
  std::copy(coords,coords+6*2,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we merge node duplicates using mergeNodes() and check values it returns.

  bool areNodesMerged; mcIdType newNbOfNodes;
  MCAuto<DataArrayIdType> arr=
    mesh->mergeNodes(0.004,areNodesMerged,newNbOfNodes);
  const mcIdType idsExpected[6] = {0, 1, 0, 2, 2, 0};
  CPPUNIT_ASSERT(std::equal(idsExpected,idsExpected+6,arr->getPointer()));
  CPPUNIT_ASSERT( areNodesMerged );
  CPPUNIT_ASSERT_EQUAL( ToIdType( 3 ), newNbOfNodes );

Contents of arr shows ids of old nodes after the merging. The nodes considered equal one to the other have the same id in arr.

Next we compare behavior of mergeNodes() and that of mergeNodes2() which, in contrast to mergeNodes(), moves merged nodes to their barycenter.
We expect that mergeNodes2() moves old nodes #0, #2 and #5 to their barycenter equal to position of node #2.
First we check that mergeNodes() does not move nodes coincident with the node #2 to the position of node #2, and then we check that mergeNodes2() does move. (We check only the second (Y) component of node coordinates since the first component of these nodes is exactly same.)

  const double* baryCoords2 = coords + 2*2; // initial coordinates of node #2
  coordsArr=mesh->getCoordinatesAndOwner(); // retrieve a new shorten coord array
  CPPUNIT_ASSERT( fabs( baryCoords2[1] - coordsArr->getIJ(0,1)) > 1e-4 ); // Y of node #0 differs from that of baryCoords2
  // restore coordinates
  coordsArr->alloc(6,2);
  std::copy(coords,coords+6*2,coordsArr->getPointer());
  mesh->setCoords(coordsArr);
  // call mergeNodesCenter()
  arr = mesh->mergeNodesCenter(0.004,areNodesMerged,newNbOfNodes);
  coordsArr=mesh->getCoordinatesAndOwner(); // retrieve a new shorten coord array
  CPPUNIT_ASSERT_DOUBLES_EQUAL( baryCoords2[1], coordsArr->getIJ(0,1), 13 ); // Y of node #0 equals to that of baryCoords2

Removing cell duplicates

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells, so that

• the cell #2 has the same nodal connectivity as the cell #1 does,
• the cell #3 has the same nodal connectivity as the cell #0 does,
• the cell #4 is based on the same nodes as the cell #0 but nodes order is different.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[11]={0,3,4,1, 1,4,2, 4,1,0,3};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+0); // 0     
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4); // 1     
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4); // 2 == 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+0); // 3 == 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+7); // 4 ~~ 0
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we use zipConnectivityTraducer() to remove duplicate cells. Then we check that two cells, having exactly same nodal connectivity with other cells, have been removed.

  const mcIdType oldNbCells = mesh->getNumberOfCells();
  DataArrayIdType *arr = mesh->zipConnectivityTraducer(0);
  CPPUNIT_ASSERT_EQUAL( oldNbCells-2, mesh->getNumberOfCells() );
  const mcIdType idsExpected[5] = {0, 1, 1, 0, 2};
  CPPUNIT_ASSERT(std::equal(idsExpected,idsExpected+5,arr->getPointer()));

Contents of arr shows ids of cells after duplicates removal. If a value (cell id) equals to its index in arr, this means that the cell is not a duplicate of any cell with lower id. Else, the value gives a cell id to which this cell is equal.
Thus, the cells #0 and #1 have no preceding equal cell since arr[i] == i.
The cell #2 equals to the cell #1 (== arr[2] ).
The cell #3 equals to the cell #0 (== arr[3] ).
The cell #4 has no equal cell. This is because the cell comparison technique specified when we called zipConnectivityTraducer() was 0 ("exact"), if we had used the technique 2 ("nodal"), arr[4] would be 0.

Removing unused nodes

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);   
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10);
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14);
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we create mesh2 including all nodes but only two cells of mesh, and we use zipCoordsTraducer() to remove unused nodes from mesh2. zipCoordsTraducer() returns an array with -1 for unused nodes and new ids for used ones.

  const mcIdType cellIds[2]={1,2};
  MEDCouplingUMesh* mesh2=(MEDCouplingUMesh*)mesh->buildPartOfMySelf(cellIds,cellIds+2,true);
  DataArrayIdType *arr=mesh2->zipCoordsTraducer();
  CPPUNIT_ASSERT_EQUAL( ToIdType(4), mesh2->getNumberOfNodes() ); // nb of nodes decreased
  CPPUNIT_ASSERT_EQUAL( mesh->getNumberOfNodes(), arr->getNumberOfTuples() );
  const mcIdType idsExpected[9] = {-1,0,1,-1,2,3,-1,-1,-1}; // -1 for unused nodes
  CPPUNIT_ASSERT(std::equal(idsExpected,idsExpected+9,arr->getPointer()));

Retrieving unused nodes

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);   
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10);
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14);
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we create mesh2 including all nodes but only two cells of mesh, and we use getNodeIdsInUse() to get nodes of mesh2 used in its two cells. getNodeIdsInUse() returns an array with -1 for unused nodes and new ids for used ones.

  const mcIdType cellIds[2]={1,2};
  MEDCouplingUMesh* mesh2=(MEDCouplingUMesh*)mesh->buildPartOfMySelf(cellIds,cellIds+2,true);
  mcIdType newNbOfNodes = 0;
  DataArrayIdType *arr=mesh2->getNodeIdsInUse( newNbOfNodes );
  const mcIdType idsExpected[9] = {-1,0,1,-1,2,3,-1,-1,-1};
  CPPUNIT_ASSERT(std::equal(idsExpected,idsExpected+9,arr->getPointer()));

Now we use newNbOfNodes returned by getNodeIdsInUse() to convert arr to "New to Old" mode.

  DataArrayIdType *arr2=arr->invertArrayO2N2N2O(newNbOfNodes);
  const mcIdType idsExpected2[4] = {1,2,4,5};
  CPPUNIT_ASSERT(std::equal(idsExpected2,idsExpected2+4,arr2->getPointer()));

Conversion of cells to "poly" types

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we convert cells #1 and #3 to type POLYGON and check the result

  const mcIdType cells[2]={1,3};
  mesh->convertToPolyTypes(cells, cells+2);
  CPPUNIT_ASSERT( mesh->getTypeOfCell(0) == INTERP_KERNEL::NORM_QUAD4 );
  CPPUNIT_ASSERT( mesh->getTypeOfCell(1) == INTERP_KERNEL::NORM_POLYGON );
  CPPUNIT_ASSERT( mesh->getTypeOfCell(2) == INTERP_KERNEL::NORM_TRI3 );
  CPPUNIT_ASSERT( mesh->getTypeOfCell(3) == INTERP_KERNEL::NORM_POLYGON );

Scaling the mesh

First, we create a 2D mesh with 4 nodes and no cells.

  double coords[4*2]={0.0,0.0, 1.0,0.0, 1.0,1.0, 0.0,1.0}; // 2D coordinates of 4 nodes
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(coords, 4,2);
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);
  DataArrayDouble *initCoords = coordsArr->deepCopy();

Then we scale it by a factor of 2 with a center (0.,0.).

  const double center[2] = {0.,0.};
  const double factor = 2.;
  mesh->scale( center, factor );

Finally we check that all node coordinates have changed by more than 0.9.

  const DataArrayDouble * coordsArr2 = mesh->getCoords();
  CPPUNIT_ASSERT( coordsArr2->isEqualWithoutConsideringStr( *initCoords, 1.0 ));
  CPPUNIT_ASSERT( !coordsArr2->isEqualWithoutConsideringStr( *initCoords, 0.9 ));
  // release data
  initCoords->decrRef();

Translating the mesh

First, we create a 2D mesh with 4 nodes and no cells.

  double coords[4*2]={0.0,0.0, 1.0,0.0, 1.0,1.0, 0.0,1.0}; // 2D coordinates of 4 nodes
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(coords, 4,2);
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);
  DataArrayDouble *initCoords = coordsArr->deepCopy();

Then we translate it by a vector (1.,1.).

  double vector[2] = {1.,1.};
  mesh->translate( vector );

Finally we check that all node coordinates have changed by more than 0.9.

  const DataArrayDouble * coordsArr2 = mesh->getCoords();
  CPPUNIT_ASSERT( coordsArr2->isEqualWithoutConsideringStr( *initCoords, 1.0 ));
  CPPUNIT_ASSERT( !coordsArr2->isEqualWithoutConsideringStr( *initCoords, 0.9 ));
  // release data
  initCoords->decrRef();

Rotating the mesh

First, we create a 2D mesh with 4 nodes and no cells.

  double coords[4*2]={0.0,0.0, 0.1,0.0, 0.1,0.1, 0.0,0.1}; // 2D coordinates of 4 nodes
  double coordsOrig[4*2];
  std::copy(coords,coords+sizeof(coords)/sizeof(double),coordsOrig);//keep tracks of initial values
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(coords, 4,2);
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);

Then we rotate it around a point (0.,0.) by 90 degrees clockwise.

  double center[3] = {0.,0.,0.}; // it suits for 2D as well
  double vector[3] = {0.,0.,1.}; // it is not used in 2D
  mesh->rotate( center, vector, -M_PI/2); // warning here C++ 'coords' array (defined above) has been modified !

Next, we make a 3D mesh from the 2D one and rotate it around the Z axis by 90 degrees counter-clockwise.

  mesh->changeSpaceDimension(3);
  mesh->rotate( center, vector, +M_PI/2);

Finally we transform the mesh back to 2D space and check that all nodes get back to the initial location.

  mesh->changeSpaceDimension(2);
  const DataArrayDouble * coordsArr2 = mesh->getCoords();
  coordsArr->useExternalArrayWithRWAccess(coordsOrig, 4,2);
  CPPUNIT_ASSERT( coordsArr2->isEqualWithoutConsideringStr( *coordsArr, 1e-13 ));

Access

Getting node coordinates along an axis

We create an 1D Cartesian mesh and retrieves node coordinates using getCoordsAt().

  const double coords[3] = {1.,2.,4.};
  MCAuto<DataArrayDouble> x = DataArrayDouble::New();
  x->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh=MEDCouplingCMesh::New();
  mesh->setCoordsAt(0,x);
  const DataArrayDouble* x2=mesh->getCoordsAt(0);
  CPPUNIT_ASSERT( x2->isEqual( *x, 1e-13 ));

Getting coordinates of a node

The following code creates a 2D MEDCouplingUMesh with 3 nodes and no cells.

  double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3};
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(coords, 3,2);
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);

Here we get coordinates of the second node and check its two coordinates.

  std::vector<double> coords2;
  mesh->getCoordinatesOfNode(1,coords2);
  CPPUNIT_ASSERT_DOUBLES_EQUAL(coords[2],coords2[0],1e-13);
  CPPUNIT_ASSERT_DOUBLES_EQUAL(coords[3],coords2[1],1e-13);

Cells correspondence in two meshes

First, we create a 2D mesh1 with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh1=MEDCouplingUMesh::New();
  mesh1->setMeshDimension(2);
  mesh1->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh1->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // #0
  mesh1->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // #1
  mesh1->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // #2
  mesh1->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // #3
  mesh1->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // #4
  mesh1->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh1->setCoords(coordsArr);

Then we create a mesh2 which includes cells #4, #2 and #0 of mesh1. The two meshes share the same node coordinates array.

  const mcIdType cells2[3] = { 4,2,0 }; // even cells selected
  MCAuto<MEDCouplingUMesh> mesh2 =
    (MEDCouplingUMesh*) mesh1->buildPartOfMySelf( cells2, cells2+3, true );

Now we ascertain that

• areCellsIncludedIn() detects that all cells of mesh2 are present in mesh1,
• the correspondence array corr2to1, which gives cell ids of mesh2 within mesh1, is equal to the array cells2 which selected cells from mesh1 for creation of mesh2.

  int compType = 0; // the strongest policy
  DataArrayIdType *corr2to1, *corr1to2;
  // a larger mesh1 includes a smaller mesh2
  CPPUNIT_ASSERT( mesh1->areCellsIncludedIn( mesh2, compType, corr2to1 ));
  CPPUNIT_ASSERT( std::equal( cells2, cells2+3, corr2to1->getConstPointer() ));

Now we apply areCellsIncludedIn() in a reverse direction and ascertain that it returns false.

  // the smaller mesh2 does NOT include the larger mesh1
  CPPUNIT_ASSERT( ! mesh2->areCellsIncludedIn( mesh1, compType, corr1to2 ));
  const mcIdType corr1to2Expected[5] = {2, 3, 1, 4, 0};
  CPPUNIT_ASSERT(std::equal( corr1to2Expected, corr1to2Expected+5, corr1to2->getConstPointer() ));

The contents of the correspondence array corr1to2 [2, 3, 1, 4, 0] means the following.

• The cell #0 of mesh1 is equal to the cell #2 (== corr1to2[ 0 ]) of mesh2.
• The cell #1 of mesh1 is missing from mesh2 (as corr1to2[ 1 ] >= mesh2->getNumberOfCells()).
• The cell #2 of mesh1 is equal to the cell #1 (== corr1to2[ 2 ]) of mesh2.
• The cell #3 of mesh1 is missing from mesh2 (as corr1to2[ 3 ] >= mesh2->getNumberOfCells()).
• The cell #4 of mesh1 is equal to the cell #0 (== corr1to2[ 4 ]) of mesh2.

Deep comparison of meshes

First, we create two 2D meshes with two triangles, so that

• their nodes are almost same but permuted,
• the first triangle is based exactly on the same nodes (taking the permutation into account),
• an order of nodes in the second triangle is changed.

  // mesh 1
  MEDCouplingUMesh *mesh1=MEDCouplingUMesh::New();
  const double coords[4*2]={0.0,0.0,  // #0
                            1.0,0.0,  // #1
                            1.0,1.0,  // #2
                            0.0,1.0}; // #3
  {
    mesh1->setMeshDimension(2);
    MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
    coordsArr->useExternalArrayWithRWAccess( coords, 4, 2 );
    mesh1->setCoords(coordsArr);
    mesh1->allocateCells(2);
    const mcIdType conn[6]={0,1,2, 1,2,3};
    mesh1->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+0);  // #0
    mesh1->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+3);  // #1
    mesh1->finishInsertingCells();
  }
  // mesh 2
  MEDCouplingUMesh *mesh2=MEDCouplingUMesh::New();
  const double coords2[4*2]={0.0,1.0,    // #0 = #3
                             0.0,0.0,    // #1 = #0
                             1.0,0.0,    // #2 = #1
                             1.0,1.001}; // #3 ~ #2
  {
    mesh2->setMeshDimension(2);
    MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
    coordsArr->useExternalArrayWithRWAccess( coords2, 4, 2 );
    mesh2->setCoords(coordsArr);
    mesh2->allocateCells(2);
    const mcIdType conn[6]={2,3,0, 3,1,2};
    mesh2->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+0);  // #0 = #1
    mesh2->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+3);  // #1 ~ #0
    mesh2->finishInsertingCells();
  }

Then we check that

• checkDeepEquivalWith() considers the meshes equal (i.e. it does not throw any exception) if it is called with a cell comparison policy cellCompPol == 1
• mapping from mesh1 to mesh2 for both nodes and cells is as expected.

  int cellCompPol = 1; // "permuted same orientation" - policy of medium severity
  DataArrayIdType *nOld2New, *cOld2New;
  mesh1->checkDeepEquivalWith( mesh2, cellCompPol, 0.002, cOld2New, nOld2New );
  const mcIdType nOld2NewExpected[4] = { 3, 0, 1, 2 };
  const mcIdType cOld2NewExpected[2] = { 1, 0 };
  CPPUNIT_ASSERT(std::equal(nOld2NewExpected,nOld2NewExpected+4,nOld2New->getConstPointer()));
  CPPUNIT_ASSERT(std::equal(cOld2NewExpected,cOld2NewExpected+2,cOld2New->getConstPointer()));

Next we ascertain that checkDeepEquivalOnSameNodesWith() consider mesh1 and mesh2 different as they do not share the same nodal connectivity array.
After that we make the meshes share the node coordinates array and insert new triangles based on the same nodes but in different order. This is to ascertain that checkDeepEquivalOnSameNodesWith() called with the weakest cell comparison policy considers the meshes equal.

  cOld2New->decrRef(); // else memory leaks
  CPPUNIT_ASSERT_THROW ( mesh1->checkDeepEquivalOnSameNodesWith( mesh2, cellCompPol, 0.002, cOld2New ), INTERP_KERNEL::Exception );
  mesh2->setCoords( mesh1->getCoords() ); // make meshes share the same coordinates array
  mesh2->allocateCells(2);
  const mcIdType conn[6]={1,2,3, 1,0,2};
  mesh2->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn+0); // #0 = #1
  mesh2->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn+3); // #1 ~ #0
  mesh2->finishInsertingCells();
  cellCompPol = 2; // the weakest policy
  mesh1->checkDeepEquivalOnSameNodesWith( mesh2, cellCompPol, 0, cOld2New );

Getting barycenters of cells

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,2,4, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we use getPartBarycenterAndOwner() to get barycenters of all but the first cell.

  const mcIdType cellIds[4] = {1,2,3,4}; // cell #0 is omitted
  MCAuto<DataArrayDouble> baryCenters=
    mesh->getPartBarycenterAndOwner( cellIds, cellIds+4 );
  CPPUNIT_ASSERT( baryCenters->getNumberOfTuples() == 4 );
  CPPUNIT_ASSERT( (int)baryCenters->getNumberOfComponents() == mesh->getSpaceDimension() );

The returned array contains 4 tuples per 2 components.

Finding cells containing a point (multi-point case)

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);   
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10);
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14);
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Then we use getCellsContainingPoints() to get cells in contact with tree points. Two of them are in contact with some cells and one is not.

  const double pos[3*2] = { 10., 10,               // point out of the mesh
                            0.3, 0.3,              // point located somewhere inside the mesh
                            coords[2], coords[3]}; // point at the node #1
  const double eps = 1e-4; // ball radius
  MCAuto<DataArrayIdType> cells, cellsIndex;
  mesh->getCellsContainingPoints( pos, 3, eps, cells, cellsIndex );
  const mcIdType cellsExpected[3]={4, 0, 1};
  const mcIdType cellsIndexExpected[4]={0, 0, 1, 3};
  CPPUNIT_ASSERT(std::equal( cellsExpected,      cellsExpected+3,      cells->begin()));
  CPPUNIT_ASSERT(std::equal( cellsIndexExpected, cellsIndexExpected+4, cellsIndex->begin()));

The contents of the result arrays cells ([4, 0, 1]) and cellsIndex ([0, 0, 1, 3]) mean the following.

• Point #0 is in contact with none (== cellsIndx[1] - cellsIndx[0]) cell.
• Point #1 is in contact with 1 (== cellsIndx[2] - cellsIndx[1]) cell whose id is #4 (== cells[ cellsIndx[ 1 ]]).
• Point #2 is in contact with 2 (== cellsIndx[3] - cellsIndx[2]) cells whose ids are #0 (== cells[ cellsIndx[ 2 ]]) and #1 (== cells[ cellsIndx[ 2 ] + 1 ]).

Finding cells containing a point

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Then we use getCellsContainingPoint() to get cells in contact with a small ball (point with precision) located near the node #4 and shifted from this node by its radius eps.

  const double* coords4  = coords + 4*2; // coordinates of the node #4
  const double eps = 1e-4; // ball radius
  const double pos[2] = { coords4[0] + eps, coords4[1] - eps }; // ball center
  std::vector<mcIdType> cellIds;
  mesh->getCellsContainingPoint( pos, eps, cellIds );
  CPPUNIT_ASSERT ( ToIdType(cellIds.size()) == mesh->getNumberOfCells() );

Since the node #4 is shared by all cells, size of the vector cellIds must be equal to the number of cells in mesh.

Getting normals of cells

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells. Orientation of the cell #1 is reversed.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we use buildPartOrthogonalField() to get normal vectors to the cells.

  const mcIdType part[4] = {1,2,3,4}; // cell #0 is omitted
  MCAuto<MEDCouplingFieldDouble> vecField=
    mesh->buildPartOrthogonalField( part, part+4 );
  CPPUNIT_ASSERT ( vecField->getArray()->getNumberOfTuples() == 4 );
  CPPUNIT_ASSERT ( vecField->getArray()->getNumberOfComponents() == 3 );

Getting volumes of cells

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells. Orientation of the cell #1 is reversed.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,2,4, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we use getPartMeasureField() to get volumes of all but the first cell. If we call getPartMeasureField() with isAbs == true, the area of the cell #1 is returned positive, else, negative that reflects its inverse orientation.

  const bool isAbs = true;
  const mcIdType part[4] = {1,2,3,4}; // cell #0 is omitted
  MCAuto<DataArrayDouble> areaArr=
    mesh->getPartMeasureField( isAbs, part, part+4 );
  CPPUNIT_ASSERT( areaArr->getIJ(0,0) > 0 ); // orientation ignored
  areaArr=mesh->getPartMeasureField( !isAbs, part, part+4 );
  CPPUNIT_ASSERT( areaArr->getIJ(0,0) < 0 ); // orientation considered
  CPPUNIT_ASSERT ( areaArr->getNumberOfTuples() == 4 );

Getting cells using the bounding box

First, we create a 2D mesh with 1 TRI3 cell. Bounding box of this cell is [0.,0., 1.,1].

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(1);
  const double coords[3*2]={0.,0., 0.,1., 1.,1};
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(coords, 3,2);
  mesh->setCoords(coordsArr);
  mesh->allocateCells(1);
  const mcIdType conn[3]={0,1,2};
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3,conn);
  mesh->finishInsertingCells();

Now we check how getCellsInBoundingBox() searches for cells using the bounding box. We use a bounding box touching the bounding box of the sole cell at one point (1.,1.).

  const double bbox[] = {1., 1., 1.001,1.001}; // xMin, xMax, yMin, yMax
  MCAuto<DataArrayIdType> cellIdsArr =
    mesh->getCellsInBoundingBox( bbox, 0.0 );
  CPPUNIT_ASSERT( cellIdsArr->getNumberOfTuples() == 0 );
  cellIdsArr = mesh->getCellsInBoundingBox( bbox, 0.1 );
  CPPUNIT_ASSERT( cellIdsArr->getNumberOfTuples() == 1 );

If getCellsInBoundingBox() is called with parameter eps == 0.0, the cell is not found because the two bounding boxes (one of the cell and the one passed as parameter) do not overlap.
If getCellsInBoundingBox() is called with parameter eps == 0.1, the cell is found because eps is used to increase the bounding box of the cell and thus the two bounding boxes intersect each other.

Getting boundary nodes

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);   
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10);
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14);
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we use findBoundaryNodes() to get ids of boundary nodes.

  MCAuto<DataArrayIdType> nodeIdsArr=mesh->findBoundaryNodes();
  CPPUNIT_ASSERT( nodeIdsArr->getNumberOfTuples() == mesh->getNumberOfNodes() - 1 );

findBoundaryNodes() returns all node ids except the node #4 which is in the middle of mesh.

Getting cells by nodes

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

In the following code we retrieve nodes of the cell #0 an then we call getCellIdsLyingOnNodes() twice with these nodes and with varying last parameter allNodes as input.

  std::vector<mcIdType> nodes;
  mesh->getNodeIdsOfCell( 0, nodes );
  const bool allNodes = true;
  DataArrayIdType* cellIdsArr1 = mesh->getCellIdsLyingOnNodes( &nodes[0], &nodes[0]+nodes.size(), allNodes);
  DataArrayIdType* cellIdsArr2 = mesh->getCellIdsLyingOnNodes( &nodes[0], &nodes[0]+nodes.size(),!allNodes);
  CPPUNIT_ASSERT_EQUAL( cellIdsArr1->getNumberOfTuples(), ToIdType( 1 ));
  CPPUNIT_ASSERT_EQUAL( cellIdsArr2->getNumberOfTuples(), mesh->getNumberOfCells() );

If the last parameter is true getCellIdsLyingOnNodes() looks for cells whose all nodes are given to it, hence it finds the cell #0 only.
If the last parameter is false getCellIdsLyingOnNodes() looks for any cell whose nodes are given to it, hence it finds all cells of mesh because all cells share the node #4.

Getting cells by nodes

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);   
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7); 
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10);
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14);
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

In the following code we retrieve nodes of two cells an then we use getCellIdsFullyIncludedInNodeIds() to find these cells by their nodes.

  const mcIdType cellIds[2]={1,2};
  std::vector<mcIdType> nodes;
  mesh->getNodeIdsOfCell( cellIds[0], nodes );
  mesh->getNodeIdsOfCell( cellIds[1], nodes );
  DataArrayIdType* cellIdsArr = mesh->getCellIdsFullyIncludedInNodeIds( &nodes[0], &nodes[0]+nodes.size());
  CPPUNIT_ASSERT(std::equal( cellIds, cellIds+2, cellIdsArr->getPointer() ));

Retrieving the descending connectivity with orientation

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we get and check the descending connectivity.

  DataArrayIdType *desc       =DataArrayIdType::New();
  DataArrayIdType *descIndx   =DataArrayIdType::New();
  DataArrayIdType *revDesc    =DataArrayIdType::New();
  DataArrayIdType *revDescIndx=DataArrayIdType::New();
  MEDCouplingUMesh * mesh2 = mesh->buildDescendingConnectivity2(desc,descIndx,revDesc,revDescIndx);
  const mcIdType descExpected[]        = {1,2,3,4,-3,5,6,7,8,-5,9,10,-2,11,12,13,-7,-10};
  const mcIdType descIndxExpected[]    = {0,4,7,10,14,18};
  const mcIdType revDescExpected[]     = {0, 0,3, 0,1, 0, 1,2, 1, 2,4, 2, 3, 3,4, 3, 4, 4};
  const mcIdType revDescIndxExpected[] = {0,1,3,5,6,8,9,11,12,13,15,16,17,18};
  CPPUNIT_ASSERT(std::equal(descExpected,descExpected+18,desc->getPointer()));
  CPPUNIT_ASSERT(std::equal(descIndxExpected,descIndxExpected+6,descIndx->getPointer()));
  CPPUNIT_ASSERT(std::equal(revDescExpected,revDescExpected+18,revDesc->getPointer()));
  CPPUNIT_ASSERT(std::equal(revDescIndxExpected,revDescIndxExpected+14,revDescIndx->getPointer()));

Here we get connectivity of the cell #2 (#3 in FORTRAN mode) of mesh2 to see how mutual orientation of cells in mesh and mesh2 is defined.

  const mcIdType cell2ConnExpect[] = {4,1};
  std::vector<mcIdType> cell2Conn;
  mesh2->getNodeIdsOfCell( 3-1, cell2Conn ); // cell #3 in FORTRAN mode
  CPPUNIT_ASSERT(std::equal(cell2ConnExpect,cell2ConnExpect+2,&cell2Conn[0]));

The contents of the result arrays desc and descIndx mean the following.

• The cell #0 of mesh (QUAD4) is bound by 4 (== descIndx[1] - descIndx[0]) segments (SEG2) of mesh2 whose ids in FORTRAN mode are
  • #1 (== desc[ descIndx[ 0 ]]),
  • #2 (== desc[ descIndx[ 0 ] + 1 ]),
  • #3 (== desc[ descIndx[ 0 ] + 2 ]) and
  • #4 (== desc[ descIndx[ 0 ] + 3 ]).
    Ids are positive since order of nodes in the corresponding cells of mesh and mesh2 are same. For example nodes of SEG2 #3 are [4,1] and nodes of QUAD4 #0 are [0,3,4,1].
• The cell #1 of mesh (TRI3) is bound by 3 (== descIndx[2] - descIndx[1]) segments of mesh2 whose ids in FORTRAN mode are:
  • #-3 (== desc[ descIndx[ 1 ]]),
  • #5 (== desc[ descIndx[ 1 ] + 1 ]) and
  • #6 (== desc[ descIndx[ 1 ] + 2 ]).
    The id -3 means that order of nodes in SEG2 #3 ([4,1]) is different from the order of these nodes in TRI3 #1: [1,4,2].
• etc.

The contents of the result arrays revDesc and revDescIndx mean the following.

• The cell #0 of mesh2 (SEG2) bounds 1 (== revDescIndx[1] - revDescIndx[0]) cell of mesh whose id is:
  • # 0 (== revDesc[ revDescIndx[ 0 ]]).
• The cell #1 of mesh2 bounds 2 (== revDescIndx[2] - revDescIndx[1]) cells of mesh whose ids are:
  • # 0 (== revDesc[ revDescIndx[ 1 ]]) and
  • # 1 (== revDesc[ revDescIndx[ 1 ] + 1 ]).
• etc.

Retrieving the descending connectivity

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we get and check the descending connectivity.

  DataArrayIdType *desc       =DataArrayIdType::New();
  DataArrayIdType *descIndx   =DataArrayIdType::New();
  DataArrayIdType *revDesc    =DataArrayIdType::New();
  DataArrayIdType *revDescIndx=DataArrayIdType::New();
  MEDCouplingUMesh * mesh2 = mesh->buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx);
  const mcIdType descExpected[]        = {0,1,2,3, 2,4,5, 6,7,4, 8,9,1,10, 11,12,6,9};
  const mcIdType descIndxExpected[]    = {0,4,7,10,14,18};
  const mcIdType revDescExpected[]     = {0, 0,3, 0,1, 0, 1,2, 1, 2,4, 2, 3, 3,4, 3, 4, 4};
  const mcIdType revDescIndxExpected[] = {0,1,3,5,6,8,9,11,12,13,15,16,17,18};
  CPPUNIT_ASSERT(std::equal(descExpected,descExpected+18,desc->getPointer()));
  CPPUNIT_ASSERT(std::equal(descIndxExpected,descIndxExpected+6,descIndx->getPointer()));
  CPPUNIT_ASSERT(std::equal(revDescExpected,revDescExpected+18,revDesc->getPointer()));
  CPPUNIT_ASSERT(std::equal(revDescIndxExpected,revDescIndxExpected+14,revDescIndx->getPointer()));

The contents of the result arrays desc and descIndx mean the following.

• The cell #0 of mesh (QUAD4) is bound by 4 (== descIndx[1] - descIndx[0]) segments (SEG2) of mesh2 whose ids are
  • #0 (== desc[ descIndx[ 0 ]]),
  • #1 (== desc[ descIndx[ 0 ] + 1 ]),
  • #2 (== desc[ descIndx[ 0 ] + 2 ]) and
  • #3 (== desc[ descIndx[ 0 ] + 3 ]).
• The cell #1 of mesh (TRI3) is bound by 3 (== descIndx[2] - descIndx[1]) segments of mesh2 whose ids are:
  • #2 (== desc[ descIndx[ 1 ]]),
  • #4 (== desc[ descIndx[ 1 ] + 1 ]) and
  • #5 (== desc[ descIndx[ 1 ] + 2 ]).
• etc.

The contents of the result arrays revDesc and revDescIndx mean the following.

• The cell #0 of mesh2 (SEG2) bounds 1 (== revDescIndx[1] - revDescIndx[0]) cell of mesh whose id is:
  • # 0 (== revDesc[ revDescIndx[ 0 ]]).
• The cell #1 of mesh2 bounds 2 (== revDescIndx[2] - revDescIndx[1]) cells of mesh whose ids are:
  • # 0 (== revDesc[ revDescIndx[ 1 ]]) and
  • # 1 (== revDesc[ revDescIndx[ 1 ] + 1 ]).
• etc.

Getting the reverse nodal connectivity

First, we create a 2D mesh with 3 QUAD4 and 2 TRI3 cells.

  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setMeshDimension(2);
  mesh->allocateCells(5);
  const mcIdType conn[18]={0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4};
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn);    // 0
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+4);  // 1
  mesh->insertNextCell(INTERP_KERNEL::NORM_TRI3,3, conn+7);  // 2
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+10); // 3
  mesh->insertNextCell(INTERP_KERNEL::NORM_QUAD4,4,conn+14); // 4
  mesh->finishInsertingCells();
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->alloc(9,2);
  const double coords[18]={-0.3,-0.3, 0.2,-0.3, 0.7,-0.3, -0.3,0.2, 0.2,0.2, 0.7,0.2, -0.3,0.7, 0.2,0.7, 0.7,0.7 };
  std::copy(coords,coords+18,coordsArr->getPointer());
  mesh->setCoords(coordsArr);

Now we get and check its reverse nodal connectivity.

  DataArrayIdType *revNodal=DataArrayIdType::New();
  DataArrayIdType *revNodalIndx=DataArrayIdType::New();
  mesh->getReverseNodalConnectivity(revNodal,revNodalIndx);
  const mcIdType revNodalExpected[18]={0,0,1,1,2,0,3,0,1,2,3,4,2,4,3,3,4,4};
  const mcIdType revNodalIndexExpected[10]={0,1,3,5,7,12,14,15,17,18};
  CPPUNIT_ASSERT(std::equal(revNodalExpected,revNodalExpected+18,revNodal->getPointer()));
  CPPUNIT_ASSERT(std::equal(revNodalIndexExpected,revNodalIndexExpected+10,revNodalIndx->getPointer()));

The contents of the result arrays mean the following.

• Node #0 is shared by 1 (== revNodalIndx[1] - revNodalIndx[0]) cell whose id is #0 (== revNodal[ revNodalIndx[ 0 ]]).
• Node #1 is shared by 2 (== revNodalIndx[2] - revNodalIndx[1]) cells whose ids are #0 (== revNodal[ revNodalIndx[ 1 ]]) and #1 (== revNodal[ revNodalIndx[ 1 ] + 1 ]).
• etc.

Getting a minimum box bounding nodes

First, we create a 3D mesh with 2 nodes, so that the first one has minimal coordinates and the second one has maximal coordinates.

  double cc[2*3]={0.0, 0.1, 0.2, // 3D coordinates of 2 nodes
                  2.0, 2.1, 2.2};
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(cc, 2,3);
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);

Now we get a bounding box enclosing these nodes. This bounding box should contain coordinates of our two nodes (but in "no interlace" mode), as the nodes coincide with points returned by the bounding box.

  double bbox[3][2];
  mesh->getBoundingBox( (double*) bbox );
 
  // check the returned coordinates of extremum points of the bounding box
  for ( mcIdType i = 0; i < 2; ++i )   // point id
    for ( mcIdType j = 0; j < 3; ++j ) // component
      CPPUNIT_ASSERT_DOUBLES_EQUAL( cc[ i*3 + j ], bbox[j][i], 1e-13);

Getting nodes close to a point

The following code creates a 2D MEDCouplingUMesh with 5 nodes and no cells.

  // 2D coordinates of 5 nodes
  double coords[5*2]={0.3,-0.30001, // #0
                      0.2,-0.3,   // #1
                      0.3,-0.30002, // #2
                      1.1,0.0,    // #3
                      0.3,-0.30003};// #4
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(coords, 5,2);
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);

Now we define an array of coordinates of a point close to nodes #0, #2 and #4.

Thus we expect that getNodeIdsNearPoint() that we are going to use, if called with eps = 0.003, would return ids of nodes #0, #2 and #4.

  double point [2]={0.3, -0.3}; // point close to nodes #0, #2 and #4
  DataArrayIdType *ids = mesh->getNodeIdsNearPoint(point, 1e-2);
 
  // check found ids
  const mcIdType expectedIDs[3] = {0,2,4};
  DataArrayIdType * okIDs = ids->findIdsEqualList ( expectedIDs, expectedIDs+3 );
  CPPUNIT_ASSERT_EQUAL(ToIdType(3), okIDs->getNumberOfTuples());
 
  // release data
  ids->decrRef();
  okIDs->decrRef();

Getting nodes close to some points

The following code creates a 2D MEDCouplingUMesh with 7 nodes and no cells.

  // 2D coordinates of 7 nodes
  double coords[7*2]={0.3,-0.301, // #0
                      0.2,-0.3,   // #1
                      0.3,-0.302, // #2
                      1.1,0.0,    // #3
                      1.1,0.0,    // #4
                      1.1,0.002,  // #5
                      0.3,-0.303};// #6
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(coords, 7,2);
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);

Now we define an array of coordinates of 3 points near which we want to find nodes of the mesh.

• Point #0 is at distance 0.001 from the node #1.
• Point #1 is rather far from all nodes.
• Point #2 is close to nodes #3, #4 and #5.

Thus we expect that getNodeIdsNearPoints() that we are going to use, if called with eps = 0.003, would return ids of close nodes #1, #3, #4 and #5.

  const mcIdType nbOfPoints = 3;
  double points [nbOfPoints*2]={0.2,-0.30001,  // ~ node #1
                                0.0, 0.0,
                                1.1, 0.002}; // ~ nodes #3, #4 and #5
  DataArrayIdType *ids, *idsIndex;
  mesh->getNodeIdsNearPoints(points, nbOfPoints, 1e-1,ids,idsIndex);
 
  // check found ids (i.e. contents of 'ids' array)
  const mcIdType expectedIDs[4] = {1, 3, 4, 5};
  DataArrayIdType * okIDs = ids->findIdsEqualList ( expectedIDs, expectedIDs+4 );
  CPPUNIT_ASSERT_EQUAL(ToIdType(4), okIDs->getNumberOfTuples());
 
  // release data
  ids->decrRef();
  idsIndex->decrRef();
  okIDs->decrRef();

idsIndex returns [0, 1, 1, 4] which means that:

• Point #0 is close to 1 (== idsIndex[1] - idsIndex[0]) node whose id is ids[ idsIndex[ 0 ]].
• Point #1 is close to 0 (== idsIndex[2] - idsIndex[1]) nodes.
• Point #2 is close to 3 (== idsIndex[3] - idsIndex[2]) nodes whose ids are ids[ idsIndex[ 2 ]], ids[ idsIndex[ 2 ] + 1 ] and ids[ idsIndex[ 2 ] + 2 ].

Finding coincident nodes

First, we create a mesh with 6 nodes, of which two nodes (#3 and #4) are fully coincident and 3 nodes (#0, #2 and #5) have distance less than 0.004 between them.

  double coords[6*2]={0.3,-0.301, // 0
                      0.2,-0.3,   // 1
                      0.3,-0.302, // 2
                      1.1,0.0,    // 3
                      1.1,0.0,    // 4
                      0.3,-0.303};// 5
  MCAuto<DataArrayDouble> coordsArr=DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess(coords, 6,2);
  MCAuto<MEDCouplingUMesh> mesh=MEDCouplingUMesh::New();
  mesh->setCoords(coordsArr);

Then, we use findCommonNodes() to find coincident nodes, and check that (1) calling findCommonNodes() with prec == 1e-13 finds the two fully coincident nodes only and (2) findCommonNodes(0.004) finds 5 equal nodes.

  DataArrayIdType *com, *comI;
  mesh->findCommonNodes(1e-13,-1,com,comI);
  CPPUNIT_ASSERT_EQUAL(ToIdType(2),com->getNumberOfTuples());
  com->decrRef(); comI->decrRef();
  mesh->findCommonNodes(0.004,-1,com,comI);
  CPPUNIT_ASSERT_EQUAL(ToIdType(5), com->getNumberOfTuples());

Arrays

Create

Building an array from scratch in C++

In this example we will create arrays with 12 tuples constituted each of 3 components. These arrays will be created using different ways.
The following code is only based using DataArrayDouble but the use of DataArrayInt is strictly equivalent.

  const mcIdType nbOfNodes=12;
  double coords[3*nbOfNodes]={2.,3.,4.,3.,4.,5.,4.,5.,6.,5.,6.,7.,6.,7.,8.,7.,8.,9.,8.,9.,10.,9.,10.,11.,10.,11.,12.,11.,12.,13.,12.,13.,14.,13.,14.,15.};
  //
  MEDCoupling::DataArrayDouble *coordsArr=0;
  double *tmp=0;

Building an array from scratch in C++, no copy no ownership

  coordsArr=MEDCoupling::DataArrayDouble::New();
  coordsArr->useArray(coords,false,MEDCoupling::DeallocType::CPP_DEALLOC,nbOfNodes,3);
  //now use coordsArr as you need
  //...
  //coordsArr is no more useful here : release it
  coordsArr->decrRef();

Building an array from scratch in C++, no copy with C++ ownership

  coordsArr=MEDCoupling::DataArrayDouble::New();
  tmp=new double[3*nbOfNodes];
  std::copy(coords,coords+3*nbOfNodes,tmp);
  coordsArr->useArray(tmp,true,MEDCoupling::DeallocType::CPP_DEALLOC,nbOfNodes,3);
  //now use coordsArr as you need
  //...
  //coordsArr is no more useful, release it
  coordsArr->decrRef();

Building an array from scratch in C++, no copy with C ownership

  coordsArr=MEDCoupling::DataArrayDouble::New();
  tmp=(double *)malloc(3*nbOfNodes*sizeof(double));
  std::copy(coords,coords+3*nbOfNodes,tmp);
  coordsArr->useArray(tmp,true,MEDCoupling::DeallocType::C_DEALLOC,nbOfNodes,3);
  //now use coordsArr as you need
  //...
  //coordsArr is no more useful here : release it
  coordsArr->decrRef();

Building an array from scratch in C++, with copy

  coordsArr=MEDCoupling::DataArrayDouble::New();
  coordsArr->alloc(nbOfNodes,3);
  tmp=coordsArr->getPointer();
  std::copy(coords,coords+3*nbOfNodes,tmp);
  coordsArr->declareAsNew();//you have modified data pointed by internal pointer notify object
  //now use coordsArr as you need
  //...
  //coordsArr is no more useful here : release it
  coordsArr->decrRef();

Copy DataArrays

Deep copy

To perform a deep copy of a DataArray instance simply invoke :

  coordsArrCpy=coordsArr->deepCopy();

or :

  coordsArrCpy=coordsArr->performCopyOrIncrRef(true);

coordsArrCpy is the deep copy of coordsArr so they are independent and their raw data has been deeply copied.

So it leads to the following behaviour :

  CPPUNIT_ASSERT(coordsArrCpy->isEqual(*coordsArr,1e-12));
  coordsArrCpy->setIJ(0,0,1000.);
  CPPUNIT_ASSERT(!coordsArrCpy->isEqual(*coordsArr,1e-12));//coordsArrCpy only has been modified

As coordsArrCpy is a copy object it needs to be deallocated in C++ like coordsArr.

  coordsArrCpy->decrRef();

Shallow copy

To perform a shallow copy of a DataArray instance simply invoke :

  coordsArrCpy=coordsArr->performCopyOrIncrRef(false);

coordsArrCpy is the shallow copy of coordsArr so they share the same raw data. In reality they are the same object. So it leads to the following behaviour to compare with the deep copy :

  CPPUNIT_ASSERT(coordsArrCpy->isEqual(*coordsArr,1e-12));
  coordsArrCpy->setIJ(0,0,1000.);
  CPPUNIT_ASSERT(coordsArrCpy->isEqual(*coordsArr,1e-12));//coordsArr and coordsArrCpy have been modified simultaneously

So here the content of coordsArr and coordsArrCpy are linked, contrary to the deep copy case.

As coordsArrCpy is a copy object, in C++, it needs to be deallocated.

  coordsArrCpy->decrRef();

Assignation by deep copy of DataArray

We start by building a instance of MEDCoupling::DataArrayDouble allocated or not. Here, instance is not allocated, only built empty.

  coordsArrCpy=MEDCoupling::DataArrayDouble::New();

Then, coordsArrCpy is assigned with the content of coordsArr.

  coordsArrCpy->deepCopyFrom(*coordsArr);

Then coordsArrCpy is a deep copy of coordsArr except that the instance of MEDCoupling::DataArrayDouble is those specified. But the behaviour is the same than those seen for deep copy.

  CPPUNIT_ASSERT(coordsArrCpy->isEqual(*coordsArr,1e-12));
  coordsArrCpy->setIJ(0,0,2000.);
  CPPUNIT_ASSERT(!coordsArrCpy->isEqual(*coordsArr,1e-12));//coordsArrCpy only has been modified

As always, in C++, coordsArrCpy is an object whose life cycle is fully independent from coordsArr so decrement is needed.

  coordsArrCpy->decrRef();

Building a permutation array

Here we create two arrays containing same values but in different order and then we use DataArrayInt::buildPermutationArr() to get an array showing in what places the values of b array are located in a array.

  DataArrayInt *a=DataArrayInt::New();
  const mcIdType vala[5]={4,5,6,7,8};
  a->alloc(5,1);
  std::copy(vala,vala+5,a->getPointer());
  DataArrayInt *b=DataArrayInt::New();
  const mcIdType valb[5]={5,4,8,6,7};
  b->alloc(5,1);
  std::copy(valb,valb+5,b->getPointer());
  DataArrayIdType *c=a->buildPermutationArr(*b);

The result array c contains [1,0,4,2,3].

Modify

Inverting renumbering maps

invertArrayO2N2N2O()

In this example we create a DataArrayInt containing a renumbering map in "Old to New" mode, convert it into the renumbering map in "New to Old" mode and check the result.

  const mcIdType arr1[6]={2,0,4,1,5,3};
  DataArrayInt *da=DataArrayInt::New();
  da->alloc(6,1);
  std::copy(arr1,arr1+6,da->getPointer());
  DataArrayIdType *da2=da->invertArrayO2N2N2O(6);
  const mcIdType expected1[6]={1,3,0,5,2,4};
  for(mcIdType i=0;i<6;i++)
    CPPUNIT_ASSERT_EQUAL(expected1[i],da2->getIJ(i,0));

invertArrayN2O2O2N()

In this example we create a DataArrayInt containing a renumbering map in "New to Old" mode, convert it into the renumbering map in "Old to New" mode and check the result.

  const mcIdType arr1[6]={2,0,4,1,5,3};
  DataArrayInt *da=DataArrayInt::New();
  da->alloc(6,1);
  std::copy(arr1,arr1+6,da->getPointer());
  DataArrayIdType *da2=da->invertArrayN2O2O2N(6);
  const mcIdType expected1[6]={1,3,0,5,2,4};
  for(mcIdType i=0;i<6;i++)
    CPPUNIT_ASSERT_EQUAL(expected1[i],da2->getIJ(i,0));

Set part of values of DataArrayDouble

setSelectedComponents()

First, we create a 'source' array.

        array1=[1.,2., 3.,4., 5.,6.]
        da=DataArrayDouble(array1,3,2)
        da.setInfoOnComponents( ["a1","a2"])

Now we create a larger zero array and assign the array da into it.

        dv=DataArrayDouble()
        dv.alloc( 4, 4 )
        dv.fillWithZero()
        dv.setInfoOnComponents( ["v1","v2","v3","v4"])
        dv2 = dv.deepCopy()
        dv.setSelectedComponents( da, [1,0] )

As result contents of the array dv are as follows.

Info of components : "a2"   "a1"   "v3"   "v4"
    Tuple #0 : 2 1 0 0
    Tuple #1 : 4 3 0 0
    Tuple #2 : 6 5 0 0
    Tuple #3 : 0 0 0 0

The same result can be achieved other way (except that component info is not copied):

        dv2[:3,[1,0]] = da
        self.assertTrue( dv.isEqualWithoutConsideringStr( dv2, 1e-20 ))

setPartOfValues1()

We create two arrays:

• a "large" (4x4) zero array da to assign to and
• a smaller (2x2) array dv filled with values [7.,8.,9.,10].

        da=DataArrayDouble()
        da.alloc( 4, 4 )
        da.setInfoOnComponents( ["v1","v2","v3","v4"])
        #
        dv=DataArrayDouble()
        dv.alloc( 4, 1 )
        dv.iota(7)
        dv.rearrange( 2 )
        dv.setInfoOnComponents( ["a1","a2"])

Now we copy dv to the middle of da.

        da.fillWithZero()
        da.setPartOfValues1( dv, 1,3,1, 1,3,1, True )

As result contents of the array da are as follows.

    Info of components :"v1"   "v2"   "v3"   "v4"
    Tuple #0 : 0 0 0 0
    Tuple #1 : 0 7 8 0
    Tuple #2 : 0 9 10 0
    Tuple #3 : 0 0 0 0

Here we re-fill da with zeros and copy dv into a component of da.

Note that the last parameter strictCompoCompare should be False in this case, else MEDCoupling::DataArrayDouble::setPartOfValues1() throws an exception because da has 2 components but only one target component is specified.

        da.fillWithZero()
        da.setPartOfValues1( dv, 0,4,1, 1,2,1, False )

    Tuple #0 : 0 7 0 0
    Tuple #1 : 0 8 0 0
    Tuple #2 : 0 9 0 0
    Tuple #3 : 0 10 0 0

Below more two variants of location of target values are shown.

        da.fillWithZero()
        da.setPartOfValues1( dv, 1,2,1, 0,4,1, False )

    Tuple #0 : 0 0 0 0
    Tuple #1 : 7 8 9 10
    Tuple #2 : 0 0 0 0
    Tuple #3 : 0 0 0 0

        da.fillWithZero()
        da.setPartOfValues1( dv, 0,3,2, 1,4,2, True )

    Tuple #0 : 0 7 0 8
    Tuple #1 : 0 0 0 0
    Tuple #2 : 0 9 0 10
    Tuple #3 : 0 0 0 0

The same result can be achieved other way:

        da2 = da.deepCopy()
        da2.fillWithZero()
        da2[ 0:3:2, 1:4:2 ] = dv
        self.assertTrue( da.isEqual( da2, 1e-20 ))

setPartOfValuesSimple1()

We create an array (4x4) da to assign to and define a value dv to assign.

        da=DataArrayDouble()
        da.alloc( 4, 4 )
        dv = 7

Now we assign dv to the middle of da.

        da.fillWithZero()
        da.setPartOfValuesSimple1( dv, 1,3,1, 1,3,1 )

As result contents of the array da are as follows.

    Tuple #0 : 0 0 0 0
    Tuple #1 : 0 7 7 0
    Tuple #2 : 0 7 7 0
    Tuple #3 : 0 0 0 0

Here we re-fill da with zeros and assign dv to a component of da.

        da.fillWithZero()
        da.setPartOfValuesSimple1( dv, 0,4,1, 1,2,1 )

    Tuple #0 : 0 7 0 0
    Tuple #1 : 0 7 0 0
    Tuple #2 : 0 7 0 0
    Tuple #3 : 0 7 0 0

Below more two variants of location of target values are shown.

        da.fillWithZero()
        da.setPartOfValuesSimple1( dv, 1,2,1, 0,4,1 )

    Tuple #0 : 0 0 0 0
    Tuple #1 : 7 7 7 7
    Tuple #2 : 0 0 0 0
    Tuple #3 : 0 0 0 0

        da.fillWithZero()
        da.setPartOfValuesSimple1( dv, 0,3,2, 1,4,2 )

    Tuple #0 : 0 7 0 7
    Tuple #1 : 0 0 0 0
    Tuple #2 : 0 7 0 7
    Tuple #3 : 0 0 0 0

The same result can be achieved other way:

        da2 = da.deepCopy()
        da2.fillWithZero()
        da2[ 0:3:2, 1:4:2 ] = dv
        self.assertTrue( da.isEqual( da2, 1e-20 ))

setPartOfValuesSimple2()

We create an array (4x4) da to assign to and define a value dv to assign.

        da=DataArrayDouble()
        da.alloc( 4, 4 )
        dv = 7

Now we assign dv to the middle of da. We explicitly specify tuples and component to assign to by a list [1,2].

        da.fillWithZero()
        da[[1,2], [1,2]] = dv

As result contents of the array da are as follows.

    Tuple #0 : 0 0 0 0
    Tuple #1 : 0 7 7 0
    Tuple #2 : 0 7 7 0
    Tuple #3 : 0 0 0 0

Here we re-fill da with zeros and assign dv to a component of da.

        da.fillWithZero()
        da[[0,1,2,3], [1]] = dv

    Tuple #0 : 0 7 0 0
    Tuple #1 : 0 7 0 0
    Tuple #2 : 0 7 0 0
    Tuple #3 : 0 7 0 0

Below more two variants of location of target values are shown.

        da.fillWithZero()
        da[[1], [0,1,2,3]] = dv

    Tuple #0 : 0 0 0 0
    Tuple #1 : 7 7 7 7
    Tuple #2 : 0 0 0 0
    Tuple #3 : 0 0 0 0

        da.fillWithZero()
        da[[0,2], [1,3]] = dv

    Tuple #0 : 0 7 0 7
    Tuple #1 : 0 0 0 0
    Tuple #2 : 0 7 0 7
    Tuple #3 : 0 0 0 0

Note

MEDCoupling::DataArrayDouble::setPartOfValuesSimple2() can't be explicitly called in Python.

setPartOfValuesSimple3()

We create an array (4x4) da to assign to and define a value dv to assign.

        da=DataArrayDouble()
        da.alloc( 4, 4 )
        dv = 7

Now we assign dv to the middle of da. We explicitly specify tuples to assign to by a list [1,2]. And we specify components to assign to using slicing: 1:3.

        da.fillWithZero()
        da[[1,2], 1:3] = dv

As result contents of the array da are as follows.

    Tuple #0 : 0 0 0 0
    Tuple #1 : 0 7 7 0
    Tuple #2 : 0 7 7 0
    Tuple #3 : 0 0 0 0

Here we re-fill da with zeros and assign dv to a component of da.

        da.fillWithZero()
        da[[0,1,2,3], 1:2] = dv

    Tuple #0 : 0 7 0 0
    Tuple #1 : 0 7 0 0
    Tuple #2 : 0 7 0 0
    Tuple #3 : 0 7 0 0

Below more two variants of location of target values are shown.

        da.fillWithZero()
        da[[1], 0:4] = dv

    Tuple #0 : 0 0 0 0
    Tuple #1 : 7 7 7 7
    Tuple #2 : 0 0 0 0
    Tuple #3 : 0 0 0 0

        da.fillWithZero()
        da[[0,2], 1:4:2] = dv

    Tuple #0 : 0 7 0 7
    Tuple #1 : 0 0 0 0
    Tuple #2 : 0 7 0 7
    Tuple #3 : 0 0 0 0

Note

MEDCoupling::DataArrayDouble::setPartOfValuesSimple3() can't be explicitly called in Python.

setPartOfValues2()

We create two arrays:

• a "large" (4x7) zero array da to assign to,
• a smaller (3x2) array dv filled with values [7.,8.,9.,10.,11.,12.].

        da=DataArrayDouble()
        da.alloc( 4, 7 )
        #
        dv=DataArrayDouble()
        dv.alloc( 6, 1 )
        dv.iota(7)
        dv.rearrange( 2 )

Now we assign the two components of dv to the components of da with indices [1,3], and the 3 tuples of dv to the 3 tuples of da with indices [0,1,2]. This is the first mode of usage.

        da.fillWithZero()
        da[ [0,1,2], [1,3] ] = dv

As result contents of the array da are as follows.

    Tuple #0 : 0  7  0  8  0  0  0
    Tuple #1 : 0  9  0 10  0  0  0
    Tuple #2 : 0 11  0 12  0  0  0
    Tuple #3 : 0  0  0  0  0  0  0

Every value of dv has been assigned to its own location within da.

Now we re-fill da with zeros and rearrange dv to have 6 components. And we assign dv to the tuples of da with indices [0,2,3] . This is the second mode of usage.

        da.fillWithZero()
        dv.rearrange( 6 )
        da[ [0,2,3], [0,2,3,4,5,6]] = dv

The contents of dv have been assigned to each of specified tuples of da. Every value of dv is repeated in the 3 specified tuples within da.

    Tuple #0 : 7  0  8  9 10 11 12
    Tuple #1 : 0  0  0  0  0  0  0
    Tuple #2 : 7  0  8  9 10 11 12
    Tuple #3 : 7  0  8  9 10 11 12

Note

MEDCoupling::DataArrayDouble::setPartOfValues2() can't be explicitly called in Python.

setPartOfValues3()

We create two arrays:

• a "large" (4x7) zero array da to assign to,
• a smaller (3x2) array dv filled with values [7.,8.,9.,10.,11.,12.].

        da=DataArrayDouble()
        da.alloc( 4, 7 )
        #
        dv=DataArrayDouble()
        dv.alloc( 6, 1 )
        dv.iota(7)
        dv.rearrange( 2 )

Now we assign the two components of dv to the components of da with indices [1,3], and the 3 tuples of dv to the 3 tuples of da with indices [0,1,2] which are specified using slicing: "0:3". This is the first mode of usage.

        da.fillWithZero()
        da[ 0:3, [1,3] ] = dv

As result contents of the array da are as follows.

    Tuple #0 : 0  7  0  8  0  0  0
    Tuple #1 : 0  9  0 10  0  0  0
    Tuple #2 : 0 11  0 12  0  0  0
    Tuple #3 : 0  0  0  0  0  0  0

Every value of dv has been assigned to its own location within da.

Now we re-fill da with zeros and rearrange dv to have 6 components. And we assign dv to the tuples of da with indices [0,2] using slice notation "0:4:2". This is the second mode of usage.

        da.fillWithZero()
        dv.rearrange( 6 )
        da[ 0:4:2, [0,2,3,4,5,6]] = dv

The contents of dv have been assigned to each of specified tuples of da. Every value of dv is repeated in the 3 specified tuples within da.

    Tuple #0 : 7  0  8  9 10 11 12
    Tuple #1 : 0  0  0  0  0  0  0
    Tuple #2 : 7  0  8  9 10 11 12
    Tuple #3 : 0  0  0  0  0  0  0

Note

MEDCoupling::DataArrayDouble::setPartOfValues3() can't be explicitly called in Python.

Set part of values of DataArrayInt

setSelectedComponents()

First, we create a 'source' array.

        da=DataArrayInt()
        array1=[1,2, 3,4, 5,6]
        da.setValues(array1,3,2)
        da.setInfoOnComponents( ["a1","a2"])

Now we create a larger zero array and assign the array da to it.

        dv=DataArrayInt()
        dv.alloc( 4, 4 )
        dv.fillWithZero()
        dv.setInfoOnComponents( ["v1","v2","v3","v4"])
        dv2 = dv.deepCopy()
        dv.setSelectedComponents( da, [1,0] )

As result contents of the array dv are as follows.

Info of components : "a2"   "a1"   "v3"   "v4"
    Tuple #0 : 2 1 0 0
    Tuple #1 : 4 3 0 0
    Tuple #2 : 6 5 0 0
    Tuple #3 : 0 0 0 0

The same result can be achieved other way (except that component info is not copied):

        dv2[:3,[1,0]] = da
        self.assertTrue( dv.isEqualWithoutConsideringStr( dv2 ))

setPartOfValues1()

We create two arrays:

• a "large" (4x4) zero array da to assign to, and
• a smaller (2x2) array dv filled with values [7,8,9,10].

        da=DataArrayInt()
        da.alloc( 4, 4 )
        da.setInfoOnComponents( ["v1","v2","v3","v4"])
        #
        dv=DataArrayInt()
        dv.alloc( 4, 1 )
        dv.iota(7)
        dv.rearrange( 2 )
        dv.setInfoOnComponents( ["a1","a2"])

Now we copy dv to the middle of da.

        da.fillWithZero()
        da.setPartOfValues1( dv, 1,3,1, 1,3,1, True )

As result contents of the array da are as follows.

    Info of components :"v1"   "v2"   "v3"   "v4"
    Tuple #0 : 0 0 0 0
    Tuple #1 : 0 7 8 0
    Tuple #2 : 0 9 10 0
    Tuple #3 : 0 0 0 0

Here we re-fill da with zeros and copy dv into a component of da.

Note that the last parameter strictCompoCompare should be False in this case, else MEDCoupling::DataArrayInt::setPartOfValues1() throws an exception because da has 2 components but only one target component is specified.

        da.fillWithZero()
        da.setPartOfValues1( dv, 0,4,1, 1,2,1, False )

    Tuple #0 : 0 7 0 0
    Tuple #1 : 0 8 0 0
    Tuple #2 : 0 9 0 0
    Tuple #3 : 0 10 0 0

Below more two variants of location of target values are shown.

        da.fillWithZero()
        da.setPartOfValues1( dv, 1,2,1, 0,4,1, False )

    Tuple #0 : 0 0 0 0
    Tuple #1 : 7 8 9 10
    Tuple #2 : 0 0 0 0
    Tuple #3 : 0 0 0 0

        da.fillWithZero()
        da.setPartOfValues1( dv, 0,3,2, 1,4,2, True )

    Tuple #0 : 0 7 0 8
    Tuple #1 : 0 0 0 0
    Tuple #2 : 0 9 0 10
    Tuple #3 : 0 0 0 0

The same result can be achieved other way:

        da2 = da.deepCopy()
        da2.fillWithZero()
        da2[ 0:3:2, 1:4:2 ] = dv
        self.assertTrue( da.isEqual( da2 ))

setPartOfValuesSimple1()

We create an array (4x4) da to assign to and define a value dv to assign.

        da=DataArrayInt()
        da.alloc( 4, 4 )
        dv = 7

Now we assign dv to the middle of da.

        da.fillWithZero()
        da.setPartOfValuesSimple1( dv, 1,3,1, 1,3,1 )

As result contents of the array da are as follows.

    Tuple #0 : 0 0 0 0
    Tuple #1 : 0 7 7 0
    Tuple #2 : 0 7 7 0
    Tuple #3 : 0 0 0 0

Here we re-fill da with zeros and assign dv to a component of da.

        da.fillWithZero()
        da.setPartOfValuesSimple1( dv, 0,4,1, 1,2,1 )

    Tuple #0 : 0 7 0 0
    Tuple #1 : 0 7 0 0
    Tuple #2 : 0 7 0 0
    Tuple #3 : 0 7 0 0

Below more two variants of location of target values are shown.

        da.fillWithZero()
        da.setPartOfValuesSimple1( dv, 1,2,1, 0,4,1 )

    Tuple #0 : 0 0 0 0
    Tuple #1 : 7 7 7 7
    Tuple #2 : 0 0 0 0
    Tuple #3 : 0 0 0 0

        da.fillWithZero()
        da.setPartOfValuesSimple1( dv, 0,3,2, 1,4,2 )

    Tuple #0 : 0 7 0 7
    Tuple #1 : 0 0 0 0
    Tuple #2 : 0 7 0 7
    Tuple #3 : 0 0 0 0

The same result can be achieved other way:

        da2 = da.deepCopy()
        da2.fillWithZero()
        da2[ 0:3:2, 1:4:2 ] = dv
        self.assertTrue( da.isEqual( da2 ))

setPartOfValuesSimple2()

We create an array (4x4) da to assign to and define a value dv to assign.

        da=DataArrayInt()
        da.alloc( 4, 4 )
        dv = 7

Now we assign dv to the middle of da. We explicitly specify tuples and component to assign to by a list [1,2].

        da.fillWithZero()
        da[[1,2], [1,2]] = dv

As result contents of the array da are as follows.

    Tuple #0 : 0 0 0 0
    Tuple #1 : 0 7 7 0
    Tuple #2 : 0 7 7 0
    Tuple #3 : 0 0 0 0

Here we re-fill da with zeros and assign dv to a component of da.

        da.fillWithZero()
        da[[0,1,2,3], [1]] = dv

    Tuple #0 : 0 7 0 0
    Tuple #1 : 0 7 0 0
    Tuple #2 : 0 7 0 0
    Tuple #3 : 0 7 0 0

Below more two variants of location of target values are shown.

        da.fillWithZero()
        da[[1], [0,1,2,3]] = dv

    Tuple #0 : 0 0 0 0
    Tuple #1 : 7 7 7 7
    Tuple #2 : 0 0 0 0
    Tuple #3 : 0 0 0 0

        da.fillWithZero()
        da[[0,2], [1,3]] = dv

    Tuple #0 : 0 7 0 7
    Tuple #1 : 0 0 0 0
    Tuple #2 : 0 7 0 7
    Tuple #3 : 0 0 0 0

Note

MEDCoupling::DataArrayInt::setPartOfValuesSimple2() can't be explicitly called in Python.

setPartOfValuesSimple3()

We create an array (4x4) da to assign to and define a value dv to assign.

        da=DataArrayInt()
        da.alloc( 4, 4 )
        dv = 7

Now we assign dv to the middle of da. We explicitly specify tuples to assign to by a list [1,2]. And we specify components to assign to using slicing: 1:3.

        da.fillWithZero()
        da[[1,2], 1:3] = dv

As result contents of the array da are as follows.

    Tuple #0 : 0 0 0 0
    Tuple #1 : 0 7 7 0
    Tuple #2 : 0 7 7 0
    Tuple #3 : 0 0 0 0

Here we re-fill da with zeros and assign dv to a component of da.

        da.fillWithZero()
        da[[0,1,2,3], 1:2] = dv

    Tuple #0 : 0 7 0 0
    Tuple #1 : 0 7 0 0
    Tuple #2 : 0 7 0 0
    Tuple #3 : 0 7 0 0

Below more two variants of location of target values are shown.

        da.fillWithZero()
        da[[1], 0:4] = dv

    Tuple #0 : 0 0 0 0
    Tuple #1 : 7 7 7 7
    Tuple #2 : 0 0 0 0
    Tuple #3 : 0 0 0 0

        da.fillWithZero()
        da[[0,2], 1:4:2] = dv

    Tuple #0 : 0 7 0 7
    Tuple #1 : 0 0 0 0
    Tuple #2 : 0 7 0 7
    Tuple #3 : 0 0 0 0

Note

MEDCoupling::DataArrayInt::setPartOfValuesSimple3() can't be explicitly called in Python.

setPartOfValues2()

We create two arrays:

• a "large" (4x7) zero array da to assign to,
• a smaller (3x2) array dv filled with values [7,8,9,10,11,12].

        da=DataArrayInt()
        da.alloc( 4, 7 )
        #
        dv=DataArrayInt()
        dv.alloc( 6, 1 )
        dv.iota(7)
        dv.rearrange( 2 )

Now we assign the two components of dv to the components of da with indices [1,3], and the 3 tuples of dv to the 3 tuples of da with indices [0,1,2]. This is the first mode of usage.

        da.fillWithZero()
        da[ [0,1,2], [1,3] ] = dv

As result contents of the array da are as follows.

    Tuple #0 : 0  7  0  8  0  0  0
    Tuple #1 : 0  9  0 10  0  0  0
    Tuple #2 : 0 11  0 12  0  0  0
    Tuple #3 : 0  0  0  0  0  0  0

Every value of dv has been assigned to its own location within da.

Now we re-fill da with zeros and rearrange dv to have 6 components. And we assign dv to the tuples of da with indices [0,2,3] . This is the second mode of usage.

        da.fillWithZero()
        dv.rearrange( 6 )
        da[ [0,2,3], [0,2,3,4,5,6]] = dv

The contents of dv have been assigned to each of specified tuples of da. Every value of dv is repeated in the 3 specified tuples within da.

    Tuple #0 : 7  0  8  9 10 11 12
    Tuple #1 : 0  0  0  0  0  0  0
    Tuple #2 : 7  0  8  9 10 11 12
    Tuple #3 : 7  0  8  9 10 11 12

Note

MEDCoupling::DataArrayInt::setPartOfValues2() can't be explicitly called in Python.

setPartOfValues3()

We create two arrays:

• a "large" (4x7) zero array da to assign to,
• a smaller (3x2) array dv filled with values [7,8,9,10,11,12].

        da=DataArrayInt()
        da.alloc( 4, 7 )
        #
        dv=DataArrayInt()
        dv.alloc( 6, 1 )
        dv.iota(7)
        dv.rearrange( 2 )

Now we assign the two components of dv to the components of da with indices [1,3], and the 3 tuples of dv to the 3 tuples of da with indices [0,1,2] which are specified using slicing: "0:3". This is the first mode of usage.

        da.fillWithZero()
        da[ 0:3, [1,3] ] = dv

As result contents of the array da are as follows.

    Tuple #0 : 0  7  0  8  0  0  0
    Tuple #1 : 0  9  0 10  0  0  0
    Tuple #2 : 0 11  0 12  0  0  0
    Tuple #3 : 0  0  0  0  0  0  0

Every value of dv has been assigned to its own location within da.

Now we re-fill da with zeros and rearrange dv to have 6 components. And we assign dv to the tuples of da with indices [0,2] using slice notation "0:4:2". This is the second mode of usage.

        da.fillWithZero()
        dv.rearrange( 6 )
        da[ 0:4:2, [0,2,3,4,5,6]] = dv

The contents of dv have been assigned to each of specified tuples of da. Every value of dv is repeated in the 3 specified tuples within da.

    Tuple #0 : 7  0  8  9 10 11 12
    Tuple #1 : 0  0  0  0  0  0  0
    Tuple #2 : 7  0  8  9 10 11 12
    Tuple #3 : 0  0  0  0  0  0  0

Note

MEDCoupling::DataArrayInt::setPartOfValues3() can't be explicitly called in Python.

Excluding coincident tuples from DataArrayDouble

The code below creates an array of real values and than an array of unique values, not closer one to another than 0.2, is retrieved from it.

        array1=[2.3,1.2,1.3,2.3,2.301,0.8]
        da=DataArrayDouble(array1,6,1)
        #
        dv=da.getDifferentValues(2e-1)
        expected2=[2.301,1.3,0.8]
        self.assertEqual(3,dv.getNbOfElems())
        for i in range(3):
            self.assertAlmostEqual(expected2[i],dv.getIJ(i,0),14)
            pass

Concatenating DataArrayDouble's by appending components

In this example we create two data arrays including same number of tuples and then we concatenate them using meldWith().

  const mcIdType sameNbTuples = 7;
 
  DataArrayDouble *da1=DataArrayDouble::New();
  da1->alloc(sameNbTuples,2);
  da1->fillWithValue(7.);
  da1->setInfoOnComponent(0,"c0da1");
  da1->setInfoOnComponent(1,"c1da1");
 
  DataArrayDouble *da2=DataArrayDouble::New();
  da2->alloc(sameNbTuples,1);
  da2->iota(0.);
  da2->setInfoOnComponent(0,"c0da2");
 
  da1->meldWith(da2);

Now the array da1 includes 7 tuples (as before) of 3 components each. Its components are: "c0da1","c1da1","c0da2".

Concatenating DataArrayInt's by appending components

In this example we create two data arrays including same number of tuples and then we concatenate them using meldWith().

  const mcIdType sameNbTuples = 7;
 
  DataArrayInt *da1=DataArrayInt::New();
  da1->alloc(sameNbTuples,2);
  da1->fillWithValue(7);
  da1->setInfoOnComponent(0,"c0da1");
  da1->setInfoOnComponent(1,"c1da1");
 
  DataArrayInt *da2=DataArrayInt::New();
  da2->alloc(sameNbTuples,1);
  da2->iota(0);
  da2->setInfoOnComponent(0,"c0da2");
 
  da1->meldWith(da2);

Now the array da1 includes 7 tuples (as before) of 3 components each. Its components are: "c0da1","c1da1","c0da2".

Access

Getting a tuple of DataArrayInt

In this simple example we create an array of integers arranged into 3 tuples per 2 components, and finally print the second tuple.

        dv=DataArrayInt()
        dv.alloc( 6, 1 )
        dv.iota(7)
        dv.rearrange( 2 )
        assert dv.getTuple( 1 ) == [9,10]

The output is

 [9, 10] 

Note that we can traverse all tuples in the array by simply iterating over it as the code below does.

        for tpl in dv:
            print(tpl)

Its output follows.

(7, 8)
(9, 10)
(11, 12)

Finding values in range in DataArrayDouble

In this example we create an array da containing same values as ones returned by range( 10 ). Then we get an array of indices of values of da being in range [ 2.5, 6 ].

  DataArrayDouble *da=DataArrayDouble::New();
  da->alloc(10,1);
  da->iota();
 
  DataArrayIdType* da2 = da->findIdsInRange( 2.5, 6 );

As result contents of the array da2 are as follows.

    Tuple #0 : 3
    Tuple #1 : 4
    Tuple #2 : 5
    Tuple #3 : 6

Finding coincident tuples in DataArrayDouble

Let's create an array of 6 tuples and 2 components that can be considered as coordinates of 6 points in 2D space.

  DataArrayDouble *da=DataArrayDouble::New();
  da->alloc(6,2);
  const double array2[12]={2.3,2.3, // 0
                           1.2,1.2, // 1
                           1.3,1.3, // 2
                           2.3,2.3, // 3
                           2.301,   // 4
                           2.301,   // 5
                           0.8,0.8};// 6
  std::copy(array2,array2+12,da->getPointer());

Now we find points that are not far each from other than 1e-1.

  DataArrayIdType *c=0,*cI=0;
  da->findCommonTuples(1.01e-1,-1,c,cI);
 
  const mcIdType expected3[5]={0,3,4,1,2};
  const mcIdType expected4[3]={0,3,5};
  CPPUNIT_ASSERT(std::equal(expected3,expected3+5,c->getConstPointer()));
  CPPUNIT_ASSERT(std::equal(expected4,expected4+3,cI->getConstPointer()));
  c->decrRef();
  cI->decrRef();
  da->decrRef();

As we can realize from the above code, a hardcoded array expected3 is equal to the raw data of a DataArrayInt c and a hardcoded array expected4 is equal to the raw data of the DataArrayInt cI.

The array c contains indices of 5 coincident points. The array cI shows us boundaries of (cI->getNumberOfTuples()-1) = 2 groups of coincident points:

• The first group starts at index 0 and includes (3 - 0) = 3 points: 0,3,4.
• The second group starts at index 3 and includes (5 - 3) = 2 points: 1,2.

Creation of a sub-part of the DataArrayDouble by selecting components

        arr1=[1.,2.,3.,4.,     # tuple 0
              11.,12.,13.,14., # tuple 1
              21.,22.,23.,24., # ...
              31.,32.,33.,34.,
              41.,42.,43.,44.]
        a1=DataArrayDouble(arr1,5,4)
        a1.setInfoOnComponent(0,"a")
        a1.setInfoOnComponent(1,"b")
        a1.setInfoOnComponent(2,"c")
        a1.setInfoOnComponent(3,"d")

We created an array a1 containing 5 tuples of 4 components each (20 values). Now we are going to create an array a2 containing some components of a1.

        arr2V=[1,2,1,2,0,0]
        a2=a1.keepSelectedComponents(arr2V)

Now each tuple of a2 includes components named "b","c","b","c","a","a". Thus the result array a2 includes 30 elements (5 tuples per 6 components).

Creation of a sub-part of the DataArrayInt by selecting components

        arr1=[1,2,3,4,     # tuple 0
              11,12,13,14, # tuple 1
              21,22,23,24, #
              31,32,33,34,
              41,42,43,44]
        a1=DataArrayInt()
        a1.setValues(arr1,5,4)
        a1.setInfoOnComponent(0,"a")
        a1.setInfoOnComponent(1,"b")
        a1.setInfoOnComponent(2,"c")
        a1.setInfoOnComponent(3,"d")

We created an array a1 containing 5 tuples of 4 components each (20 values). Now we are going to create an array a2 containing some components of a1.

        arr2V=[1,2,1,2,0,0]
        a2=a1.keepSelectedComponents(arr2V)

Now each tuple of a2 includes components named "b","c","b","c","a","a". Thus the result array a2 includes 30 elements (5 tuples per 6 components).

Note that DataArrayInt::keepSelectedComponents() is called, providing the same result, by the following python code:

        a3=a1[:,arr2V ]

Other

Input/Output

Writing fields in a VTK file

In this example we

• create an 2D mesh and 3 fields on it,
• use WriteVTK() to write all the fields and the mesh to a VTK file.

  // mesh1
  const double coords[3] = {0.,2.,4.};
  MCAuto<DataArrayDouble> coordsArr = DataArrayDouble::New();
  coordsArr->useExternalArrayWithRWAccess( coords, 3, 1 );
  MCAuto<MEDCouplingCMesh> mesh1 = MEDCouplingCMesh::New();
  mesh1->setCoords(coordsArr,coordsArr); // mesh becomes a 2D one
 
  // 3 fields (lying on the same mesh!)
  MCAuto<MEDCouplingFieldDouble> field1 =
    mesh1->getMeasureField( true );
  MCAuto<MEDCouplingFieldDouble> field2 =
    mesh1->buildOrthogonalField();
  MCAuto<MEDCouplingFieldDouble> field3 =
    mesh1->fillFromAnalytic( ON_CELLS, 1, "x");
  field2->setName( "Normal" ); //  name is necessary!
  field3->setName( "Barycenter" ); //  name is necessary!
 
  // WriteVTK
  const char fileName[] = "testExample_MEDCouplingFieldDouble_WriteVTK.vtk";
  std::vector<const MEDCouplingFieldDouble *> fs( 3 ); // field series
  fs[0] = field1;
  fs[1] = field2;
  fs[2] = field3;
  MEDCouplingFieldDouble::WriteVTK( fileName, fs );

Loading a mesh using advanced API

const char fileName[]=...
const char meshName[]=...
MEDFileUMesh *medmesh=MEDFileUMesh::New(fileName,meshName)
std::vector<int> nel=medmesh->getNonEmptyLevels()
if(nel.size()<1)
  throw INTERP_KERNEL::Exception("The test is not good for my file ! Expecting a multi level mesh to play with !")
MEDCouplingUMesh *m0=medmesh->getMeshAtLevel(nel[1],false)
MEDCouplingUMesh *g1=medmesh->getGroup(nel[1],"mesh2",false)
DataArrayInt *dag1=medmesh->getGroupArr(nel[1],"mesh2",false)
MEDCouplingUMesh *g1bis=m0->buildPartOfMySelf(dag1->getConstPointer(),dag1->getConstPointer()+dag1->getNbOfElems())
g1bis->setName(dag1->getName())
if(!g1->isEqual(g1bis,1e-12))
  throw INTERP_KERNEL::Exception("hmmmm :g1 and g1bis should be equal...")
//
dag1->decrRef()
g1->decrRef()
m0->decrRef()
medmesh->decrRef()

Writing a mesh using advanced API

MEDCouplingUMesh *m=... //m is a mesh with meshDim=2 spaceDim=2
MEDCouplingUMesh *m1=... //m1 is a mesh with meshDim=1 spaceDim=2 same coords than m
MEDCouplingUMesh *m2=... //m2 is a mesh with meshDim=0 spaceDim=2 same coords than m
MEDFileUMesh *mm=MEDFileUMesh::New()
mm->setName("mm")//name needed to be non empty
mm->setDescription("Description mm")
mm->setCoords(m1->getCoords())
mm->setMeshAtLevel(-1,m1,false)
mm->setMeshAtLevel(0,m,false)
mm->setMeshAtLevel(-2,m2,false)
DataArrayInt *g1=DataArrayInt::New()
g1->alloc(2,1)
g1->setName("G1")
const int val1[2]={1,3}
std::copy(val1,val1+2,g1->getPointer())
DataArrayInt *g2=DataArrayInt::New()
g2->alloc(3,1)
g2->setName("G2")
const int val2[3]={1,2,3}
std::copy(val2,val2+3,g2->getPointer())
//
std::vector<const DataArrayInt *> grps(2)
grps[0]=g1 grps[1]=g2
mm->setGroupsAtLevel(0,grps,false)
//
g2->decrRef()
g1->decrRef()
//
mm->write(2)