Skip to content

Commit

Permalink
Add FieldMeshIntegral get set mesh API
Browse files Browse the repository at this point in the history
For cmlibs#230.
Note FieldMeshIntegralSquares has FieldMeshIntegral as a base class so these new API also work for it.
  • Loading branch information
rchristie committed Jul 4, 2022
1 parent 49c7e0b commit 3e2ecd1
Show file tree
Hide file tree
Showing 6 changed files with 149 additions and 3 deletions.
3 changes: 3 additions & 0 deletions CHANGELOG.txt
Original file line number Diff line number Diff line change
@@ -1,5 +1,8 @@
CHANGE LOG: OpenCMISS-Zinc Library

v3.8.0
Add missing type-specific field API.

v3.7.0
Support optimisation at time for time-varying fields.
Add Region change notifier, for tree structure changes.
Expand Down
2 changes: 1 addition & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ if(POLICY CMP0086)
endif()

# This is the project name and shows up in IDEs
project(Zinc VERSION 3.7.0 LANGUAGES C CXX)
project(Zinc VERSION 3.8.0 LANGUAGES C CXX)

# Set this to the development release or release candidate version.
set(Zinc_DEVELOPER_VERSION ".dev0")
Expand Down
18 changes: 18 additions & 0 deletions src/api/opencmiss/zinc/fieldmeshoperators.h
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,24 @@ ZINC_C_INLINE cmzn_field_id cmzn_field_mesh_integral_base_cast(
ZINC_API int cmzn_field_mesh_integral_destroy(
cmzn_field_mesh_integral_id *mesh_integral_field_address);

/**
* Get the mesh integrated over.
* @param mesh_integral_field Handle to the mesh integral field to query.
* @return Handle to the mesh, or NULL/invalid handle if invalid field.
*/
ZINC_API cmzn_mesh_id cmzn_field_mesh_integral_get_mesh(
cmzn_field_mesh_integral_id mesh_integral_field);

/**
* Set the mesh integrated over.
* @param mesh_integral_field Handle to the mesh integral field to modify.
* @param mesh The mesh to integrate over.
* @return Result OK on success, or ERROR_ARGUMENT if invalid mesh or mesh is
* not from this region.
*/
ZINC_API int cmzn_field_mesh_integral_set_mesh(
cmzn_field_mesh_integral_id mesh_integral_field, cmzn_mesh_id mesh);

/**
* Get the numbers of quadrature points in each element axis.
* @see cmzn_field_mesh_integral_set_numbers_of_points
Expand Down
10 changes: 10 additions & 0 deletions src/api/opencmiss/zinc/fieldmeshoperators.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,16 @@ class FieldMeshIntegral : public Field
Field(reinterpret_cast<cmzn_field_id>(field_mesh_integral_id))
{ }

Mesh getMesh() const
{
return Mesh(cmzn_field_mesh_integral_get_mesh(this->getDerivedId()));
}

int setMesh(const Mesh& mesh) const
{
return cmzn_field_mesh_integral_set_mesh(this->getDerivedId(), mesh.getId());
}

int getNumbersOfPoints(int valuesCount, int *valuesOut) const
{
return cmzn_field_mesh_integral_get_numbers_of_points(getDerivedId(),
Expand Down
45 changes: 43 additions & 2 deletions src/computed_field/computed_field_mesh_operators.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -116,12 +116,31 @@ class Computed_field_mesh_integral : public Computed_field_core

char* get_command_string();

cmzn_mesh_id getMesh()
/** @return Non-accessed mesh */
cmzn_mesh_id getMesh() const
{
return mesh;
}

int getNumbersOfPoints(int valuesCount, int *values)
/** @return Result OK on success, ERROR_ARGUMENT if mesh missing or from wrong region */
int setMesh(cmzn_mesh *meshIn)
{
if ((!meshIn) || (cmzn_mesh_get_region_internal(meshIn) != this->field->getRegion()))
{
display_message(ERROR_MESSAGE, "FieldMeshIntegral setMesh: Invalid mesh");
return CMZN_ERROR_ARGUMENT;
}
if (!cmzn_mesh_match(this->mesh, meshIn))
{
cmzn_mesh *oldMesh = this->mesh;
this->mesh = cmzn_mesh_access(meshIn);
cmzn_mesh_destroy(&oldMesh);
this->field->setChanged();
}
return CMZN_OK;
}

int getNumbersOfPoints(int valuesCount, int *values) const
{
if ((0 == valuesCount) || ((0 < valuesCount) && values))
{
Expand Down Expand Up @@ -727,6 +746,28 @@ inline Computed_field_mesh_integral *Computed_field_mesh_integral_core_cast(
reinterpret_cast<Computed_field*>(mesh_integral_field)->core));
}

cmzn_mesh_id cmzn_field_mesh_integral_get_mesh(
cmzn_field_mesh_integral_id mesh_integral_field)
{
if (mesh_integral_field)
{
Computed_field_mesh_integral *mesh_integral_core = Computed_field_mesh_integral_core_cast(mesh_integral_field);
return cmzn_mesh_access(mesh_integral_core->getMesh());
}
return nullptr;
}

int cmzn_field_mesh_integral_set_mesh(
cmzn_field_mesh_integral_id mesh_integral_field, cmzn_mesh_id mesh)
{
if (mesh_integral_field)
{
Computed_field_mesh_integral *mesh_integral_core = Computed_field_mesh_integral_core_cast(mesh_integral_field);
return mesh_integral_core->setMesh(mesh);
}
return CMZN_ERROR_ARGUMENT;
}

int cmzn_field_mesh_integral_get_numbers_of_points(
cmzn_field_mesh_integral_id mesh_integral_field,
int values_count, int *values)
Expand Down
74 changes: 74 additions & 0 deletions tests/fieldmodule/fieldmeshoperators.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -64,12 +64,15 @@ TEST(ZincFieldMeshIntegral, quadrature)

FieldMeshIntegral volumeField = zinc.fm.createFieldMeshIntegral(integrandField, coordinateField, mesh3d);
EXPECT_TRUE(volumeField.isValid());
EXPECT_EQ(mesh3d, volumeField.getMesh());
EXPECT_EQ(Element::QUADRATURE_RULE_GAUSSIAN, volumeField.getElementQuadratureRule());
FieldMeshIntegral surfaceAreaField = zinc.fm.createFieldMeshIntegral(integrandField, coordinateField, exteriorMesh2d);
EXPECT_TRUE(surfaceAreaField.isValid());
EXPECT_EQ(exteriorMesh2d, surfaceAreaField.getMesh());
EXPECT_EQ(Element::QUADRATURE_RULE_GAUSSIAN, surfaceAreaField.getElementQuadratureRule());
FieldMeshIntegral lengthField = zinc.fm.createFieldMeshIntegral(integrandField, coordinateField, mesh1d);
EXPECT_TRUE(lengthField.isValid());
EXPECT_EQ(mesh1d, lengthField.getMesh());
EXPECT_EQ(Element::QUADRATURE_RULE_GAUSSIAN, lengthField.getElementQuadratureRule());

const double scaleFactors[3] = { 1.5, 0.75, 2.0 };
Expand Down Expand Up @@ -265,6 +268,71 @@ TEST(ZincFieldMeshIntegral, quadrature)
}
}

TEST(ZincFieldMeshIntegral, cube_setMesh)
{
ZincTestSetupCpp zinc;
int result;

EXPECT_EQ(RESULT_OK, result = zinc.root_region.readFile(
TestResources::getLocation(TestResources::FIELDMODULE_CUBE_RESOURCE)));

Mesh mesh3d = zinc.fm.findMeshByDimension(3);
EXPECT_EQ(1, mesh3d.getSize());
Mesh mesh2d = zinc.fm.findMeshByDimension(2);
EXPECT_EQ(6, mesh2d.getSize());
Mesh mesh1d = zinc.fm.findMeshByDimension(1);
EXPECT_EQ(12, mesh1d.getSize());

Field coordinates = zinc.fm.findFieldByName("coordinates");
EXPECT_TRUE(coordinates.isValid());

const double oneValue = 1.0;
Field one = zinc.fm.createFieldConstant(1, &oneValue);
EXPECT_TRUE(one.isValid());

FieldMeshIntegral meshIntegral = zinc.fm.createFieldMeshIntegral(one, coordinates, mesh3d);
EXPECT_TRUE(meshIntegral.isValid());
EXPECT_EQ(mesh3d, meshIntegral.getMesh());

Fieldcache fieldcache = zinc.fm.createFieldcache();
const double TOL = 1.0E-12;
double volume;
EXPECT_EQ(RESULT_OK, meshIntegral.evaluateReal(fieldcache, 1, &volume));
EXPECT_NEAR(1.0, volume, TOL);

EXPECT_EQ(RESULT_OK, meshIntegral.setMesh(mesh2d));
EXPECT_EQ(mesh2d, meshIntegral.getMesh());
double area;
EXPECT_EQ(RESULT_OK, meshIntegral.evaluateReal(fieldcache, 1, &area));
EXPECT_NEAR(6.0, area, TOL);

EXPECT_EQ(RESULT_OK, meshIntegral.setMesh(mesh1d));
EXPECT_EQ(mesh1d, meshIntegral.getMesh());
double length;
EXPECT_EQ(RESULT_OK, meshIntegral.evaluateReal(fieldcache, 1, &length));
EXPECT_NEAR(12.0, length, TOL);

FieldElementGroup elementGroup2 = zinc.fm.createFieldElementGroup(mesh2d);
EXPECT_TRUE(elementGroup2.isValid());
MeshGroup meshGroup2 = elementGroup2.getMeshGroup();
EXPECT_TRUE(meshGroup2.isValid());
EXPECT_EQ(0, meshGroup2.getSize());
EXPECT_EQ(RESULT_OK, meshIntegral.setMesh(meshGroup2));
EXPECT_EQ(meshGroup2, meshIntegral.getMesh());
EXPECT_EQ(RESULT_OK, meshIntegral.evaluateReal(fieldcache, 1, &area));
EXPECT_NEAR(0.0, area, TOL);

Element face2 = mesh2d.findElementByIdentifier(2);
EXPECT_EQ(RESULT_OK, meshGroup2.addElement(face2));
EXPECT_EQ(RESULT_OK, meshIntegral.evaluateReal(fieldcache, 1, &area));
EXPECT_NEAR(1.0, area, TOL);

Element face5 = mesh2d.findElementByIdentifier(5);
EXPECT_EQ(RESULT_OK, meshGroup2.addElement(face5));
EXPECT_EQ(RESULT_OK, meshIntegral.evaluateReal(fieldcache, 1, &area));
EXPECT_NEAR(2.0, area, TOL);
}

TEST(ZincFieldMeshIntegral, midpoint_quadrature_large_numbers)
{
ZincTestSetupCpp zinc;
Expand Down Expand Up @@ -349,8 +417,14 @@ TEST(ZincFieldMeshIntegralSquares, quadrature)

FieldMeshIntegralSquares integralSquaresField = zinc.fm.createFieldMeshIntegralSquares(integrandField, coordinateField, mesh3d);
EXPECT_TRUE(integralSquaresField.isValid());
EXPECT_EQ(mesh3d, integralSquaresField.getMesh());
EXPECT_EQ(Element::QUADRATURE_RULE_GAUSSIAN, integralSquaresField.getElementQuadratureRule());

// test can cast to MeshIntegral field and query it, as it's a base class of this field
FieldMeshIntegral integralField = integralSquaresField.castMeshIntegral();
EXPECT_TRUE(integralField.isValid());
EXPECT_EQ(mesh3d, integralField.getMesh());

const double expectedValue[] = { 3.0 + 1.0/6.0, 4.0*(3.0 + 1.0/6.0) };
double value[2];
const double tolerance = 1.0E-12;
Expand Down

0 comments on commit 3e2ecd1

Please sign in to comment.