Example a/curvature: Computing curvature values on the surface of a mesh

This example shows how to calculate simple curvature with respect to curves parameterised by xi values as well as how to calculate mean and gaussian curvature for a mesh. Curvature values can be used to shade the surface of a mesh to gain insight as to where the surface is curvier. The test mesh used in this example has only four cubic hermite elements with curvature in one xi direction but not the other.

Visualisation of curvature on a four element mesh

Simple curvature is a scalar value that gives an indication of the "curviness" of a curve at any point along the arc's path. One definition for curvature is the magnitude of a change in angle with respect to a change in arc length. If the angle changes quickly at a point then the arc is deemed to be curvy at that point. A good way of visualising curvature values is to imagine a circle that would match up with the small arc segment either side of a point of interest. The value of the curvature at this point is equal to 1/r where r is the radius of the matching circle.

There are also various more complicated curvatures, including values which give us an indication of the curviness of a surface, rather than just that of an arc in space. Mean and gaussian curvature are two common measures used for surfaces.

This example makes heavy use of computed fields to compute the formula for curvature values. It also includes a subroutine that uses computed fields to calculate the determinant of a 3x3 matrix field. Reading through the example comments will introduce you to the concept of curvature in more depth, including the formulas used.

Note that curvature values are significantly affected by scale. For example the curvature of a path walked on the the surface of the earth is very small when compared with the surface of a tennis ball - even though both are approximately spherical).

Created by Peter Bier, May 2007.

Screenshot of example a/curvature


The comfile run by this example is as follows:


# IMPORTANT information on using curvature with cubic hermite finite elements

# You will notice the formulae for curvature used in this example depend 
# heavily on 1st and 2nd derivatives. They are particularly sensitive to 
# the value of the 2nd derivative in the xi1 and xi2 directions.  
# Note that for a cubic hermite mesh, the 2nd derivatives in the xi1 and xi2
# directions are NOT continuous across elements. This means that the 
# curvature values calculated will also show discontinuities across elements.

# The formulae calculate curvature VERY nicely if the surface you are 
# describing can be perfectly described by cubic hermite basis functions. 
# Unfortunately this situation is pretty much non-existent in real world 
# geometry.

# It is important to bear in mind that this com file gives you curvature 
# values based on your finite element mesh which is only an approximation 
# of an actual physical object.

# If you examine the curvature of a finite element mesh for a cylinder, cone 
# or sphere you will notice the curvature results DO NOT agree perfectly 
# with the solution for an ideal cylinder, cone or sphere. The fact that 
# the results differ from the ideal solution demonstrates that a finite 
# element mesh cannot perfectly describe a cylinder, cone or sphere.  
# A cubic function will only ever give an approximate description of a 
# surface that is defined using sin and cos funcions.
#
#
# It can still be useful to compare the curvature value of one finite 
# element mesh with the values of another finite element mesh. Basically 
# look at a few different meshes and see if the pretty pictures are 
# useful! 

# We begin by reading in a mesh to examine its curvature.
# This example uses a very simple four element mesh.
# Modify the two lines below with if you want to read in your own mesh.
gfx read node $example/simple_element.exnode;
gfx read elements $example/simple_element.exelem;

# Add a subgroup in case we want to read in a more complicated mesh and 
# only wish to display a few elements
gfx cre egroup display_elements;
gfx mod egroup display_elements add 1..4

# For a curve r in 3d, defined in terms of a parameter t, the curvature
# can be calculated using first and second derivatives of r 
# with respect to t using the vector formula:
#
# k(t) = | r'(t) x r''(t)|
#        ________________
#             |r'(t)|^3
# 
# where k is the curvature at the point defined by t.
# Note that the curvature is a positive scalar value.

# This formula can also be used to extract a scalar curvature value at a point
# on a surface in 3D by looking at the curvature of a PARTICLIAR curve 
# through that point. Note that the curvature value we obtain will depend 
# on WHICH curve through the point we select.

# For example if we have a parameterised surface s(t,u) then by fixing u 
# at a constant value c , we obtain a formula for a curve r(t) = s(t,c)
# that lies on the surface.  By fixing t instead of u we would obtain a 
# different curve. There are an infinite number of curves that pass through
# a point on a surface, each with a potentially different curvature value.

# The surface described by a 2D finite element is parameterised in terms of
# xi1 and xi2.  This means the above idea can be used to extract curvature
# values for this surface by fixing one of the xi parameters.

# extract first order derivatives of the 3d coordinate field 
# with respect to xi1 and xi2
gfx def field coordinates_d1ds1 basis_derivative fe_field coordinates order 1 xi_indices 1;
gfx def field coordinates_d1ds2 basis_derivative fe_field coordinates order 1 xi_indices 2;

# extract second order derivatives of the 3d coordinate field 
# with respect to xi1 and xi2
gfx def field coordinates_d2ds11 basis_derivative fe_field coordinates order 2 xi_indices 1 1;
gfx def field coordinates_d2ds22 basis_derivative fe_field coordinates order 2 xi_indices 2 2;

# extract the cross derivatives, note these should be the same
gfx def field coordinates_d2ds12 basis_derivative fe_field coordinates order 2 xi_indices 1 2;
gfx def field coordinates_d2ds21 basis_derivative fe_field coordinates order 2 xi_indices 2 1;

# IMPORTANT: These derivative values are only defined on the elements.  
# They are NOT defined on the nodes. If you choose to display node points
# and attempt to label them with the above derivative fields (or fields that 
# depend on them) you will get weird behaviour  and cmgui may crash.
# If interested in viewing values, add labels to element points rather than
# node points.

# We can now use computed fields to calculate curvature values of interest

# Set up fields to calculate curvature for curves defined by fixing xi2
# and using xi1 as our parameter
gfx define field kxi1_A cross_product field coordinates_d1ds1 coordinates_d2ds11;
gfx define field kxi1_B magnitude field kxi1_A;
gfx define field kxi1_C magnitude field coordinates_d1ds1;
gfx define field kxi1_D power field kxi1_C [3];
gfx define field kxi1 divide fields kxi1_B kxi1_D;

# set up fields to calculate curvature for the curve defined by fixing xi1
# and using xi2 as our parameter
gfx define field kxi2_A cross_product field coordinates_d1ds2 coordinates_d2ds22;
gfx define field kxi2_B magnitude field kxi2_A;
gfx define field kxi2_C magnitude field coordinates_d1ds2;
gfx define field kxi2_D power field kxi2_C [3];
gfx define field kxi2 divide fields kxi2_B kxi2_D;

# Now at every point on the surface we have two curvature values kxi1 and kxi2
# which are the curvatures of the corresponding parameterized curves 
# on the surface that pass through that point.  

# We can view curvature information for a surface element by shading it with a 
# spectrum based on one of the the curvature values.  For our test case
# there is no curvature in the xi2 direction so we will shade the mesh 
# using the more interesting kxi1 field.

# create a window to display the mesh with curvature information
gfx create win;
gfx modify win 1 set perturb;

# create a spectrum for displaying information
gfx create spectrum curvature;

# increate the discretization to display the element shading better
gfx modify g_element simple_element general element_discretization "12*12*12";
gfx modify g_element display_elements general element_discretization "12*12*12";

# label some element points with curvature values
# (recall that we CANNOT label node points with curvature values as the field
#  is only defined on the elements - it is undefined at the nodes)
gfx modify g_element display_elements element_points as points glyph point label kxi1 discretization "5*2*3";

# draw a surface element shading it with the kxi1 data
gfx modify g_element display_elements lines;
gfx modify g_element display_elements surfaces data kxi1 spectrum curvature;

# auto range the spectrum 
gfx modify spectrum curvature autorange;


# The kxi1 values tell us something about how the curvature varies across 
# the surface with respect to curves defined by the xi parameters but it
# would be good to have a measure that is indepent of these parameters.

# For a surface two such measures are guassian curvature and mean curvature.

# The guassian and mean curvature can be defined in terms of the
# prinicpal curvatures K1 and K2.
#
# Recall that at any point on a surface, the value of curvature depends on 
# which of the infinite possible curves through that point you look at.
# At least one curve through the point will result in a maximum
# curvature (K1) while another will result in a minimum (K2).

# One way to visualize how this happens is to imagine a plane spanned
# by the normal to the surface at a point and a randomly chosen tangent vector.
# The intersection of the surface and this plane defines a curve in 3d 
# space which will have a particular curvature value at the point of interest. 
# As we rotate the tangent vector around, the plane also rotates and 
# the curvature value changes, increasing or decreasing.  
# One particular tangent vector will correspond to the maximum curvature 
# value (and hence defines one of the principal axes of curvature)
# The other axis is defined by a different tangent vector that corresponds to 
# the minimum curvature value. The two axes are perpendicular to each other.
#
# note: IF kxi1 and kxi2 are the PRINCIPAL CURVATURES then we can calculate 
# the mean and gaussian curvature as follows

gfx define field mean_curvature_test add fields kxi1 kxi2 scale_factors 0.5 0.5;
gfx define field gaussian_curvature_test multiply fields kxi1 kxi2;

# WARNING!
# Usually kxi1 and kxi2 do NOT coincide with the principal curvatures,
# For this contrived example kxi1 and kxi2 happen to coincide with the 
# principal curvatures which is useful for checking the accuracy of 
# the generic method.  In general this will NOT be the case.

# Also note that if you do not know the axes of principal curvature
# then it is very difficult to find K1 and K2.

# This means that in general the mean and gaussian curvature are calculated 
# using a much more complicated formula, involving the determinants of 
# matrices based on derivatives of the coordinate field of the surface 
# (The derivatives are with respect to xi1 and xi2).  These formulae will
# be introduced shortly.

# Calculating the determinant of a matrix can be done with computed
# fields using the multiply and add field operations.  
# Unfortunately we can not multiply or add more that two fields per
# operation so each term in the standard determinant formula requires 
# several fields to calculate it.

# As we will need to calculate a number of determinants I have written 
# a function which takes the name of a field of 9 values 
# (representing a 3x3 matrix) and produces a field containing the 
# corresponding determinant values.


# subroutine to calculate the determinant of a matrix
#          | A1 A2 A3 |
# det(A) = | A4 A5 A6 | 
#          | A7 A8 A9 |
#
#        = A1*A5*A9 + A2*A6*A7 + A3*A4*A8 - A3*A5*A7 - A2*A4*A9 - A1*A6*A8  
#
# inputs: name of field to calculate the determinant of
#         name of field to store determinant values in
sub det
{
	my $matrixA = shift;
	my $output = shift;
	
	print "Creating determinant for $matrixA as field $output\n";

	# calculate the 6 products for the determinant
	gfx def field $matrixA.15 multiply fields "$matrixA.1" "$matrixA.5"
	gfx def field $matrixA.159 multiply fields $matrixA.15 "$matrixA.9"

	gfx def field $matrixA.26 multiply fields "$matrixA.2" "$matrixA.6"
	gfx def field $matrixA.267 multiply fields $matrixA.26 "$matrixA.7"

	gfx def field $matrixA.34 multiply fields "$matrixA.3" "$matrixA.4"
	gfx def field $matrixA.348 multiply fields $matrixA.34 "$matrixA.8"

	gfx def field $matrixA.35 multiply fields "$matrixA.3" "$matrixA.5"
	gfx def field $matrixA.357 multiply fields $matrixA.35 "$matrixA.7"

	gfx def field $matrixA.24 multiply fields "$matrixA.2" "$matrixA.4"
	gfx def field $matrixA.249 multiply fields $matrixA.24 "$matrixA.9"

	gfx def field $matrixA.16 multiply fields "$matrixA.1" "$matrixA.6"
	gfx def field $matrixA.168 multiply fields $matrixA.16 "$matrixA.8"


	# sum them all up
	gfx def field $matrixA.sum1 add fields $matrixA.159 $matrixA.267 scale_factors 1 1
	gfx def field $matrixA.sum2 add fields $matrixA.sum1 $matrixA.348 scale_factors 1 1
	gfx def field $matrixA.sum3 add fields $matrixA.sum2 $matrixA.357 scale_factors 1 -1
	gfx def field $matrixA.sum4 add fields $matrixA.sum3 $matrixA.249 scale_factors 1 -1
	gfx def field $output add fields $matrixA.sum4 $matrixA.168 scale_factors 1 -1

}

# Armed with the determinant function we now have all the tools to
# calculate the mean and guassian curvature.
#
# To make the curvature formula easier to read I will write the partial 
# derivatives of the coordinates (x) with respect to the variables u and v 
# (where u corresponds to xi1 and v corresponds to xi2)
# using this notation the partial derivative of x with respect to xi1 is 
# written as x_u and the second order partial derivative of x 
# with respect to xi1 and xi2 is written as x_uv

# the formula for mean curvature is given by:
#
# K = det(x_uu x_u x_v) det(x_vv x_u x_v) - [det(x_uv x_u x_v)]^2
#     -----------------------------------------------------------
#                 [|x_u|^2 |x_v|^2 - (x_u . x_v)^2]^2
#
# the formula for gaussian curvature is given by:
#
# K = det(x_uu x_u x_v) |x_v|^2 - 2 det(x_uv x_u x_v)(x_u . x_v)   + 
#     ----------------------------------------------------------
#                2 [|x_u|^2 |x_v|^2 - (x_u . x_v)^2]^(3/2)
#
#             det(x_vv x_u x_v) |x_u|^2
#     ----------------------------------------
#     2 [|x_u|^2 |x_v|^2 - (x_u . x_v)^2]^(3/2)


# These formula are from taken from the chapter 
# "The Gaussian and Mean curvatures" from the book
# "Modern Differential Geometry of Curves and Surfaces with Mathematica" 
# by Gray, A

# extract first order derivatives of the 3d coordinate field 
# with respect to xi1 and xi2
gfx def field x_u basis_derivative fe_field coordinates order 1 xi_indices 1;
gfx def field x_v basis_derivative fe_field coordinates order 1 xi_indices 2;

# extract second order derivatives of the 3d coordinate field 
# with respect to xi1 and xi2
gfx def field x_uu basis_derivative fe_field coordinates order 2 xi_indices 1 1;
gfx def field x_vv basis_derivative fe_field coordinates order 2 xi_indices 2 2;

# extract cross derivatives, should be the same, so we only really 
# need the first field.
gfx def field x_uv basis_derivative fe_field coordinates order 2 xi_indices 1 2;
gfx def field x_vu basis_derivative fe_field coordinates order 2 xi_indices 2 1;

# set up matrix fields for gaussian curvature calculation
gfx define field matrixA composite x_uu x_u x_v;
gfx define field matrixB composite x_vv x_u x_v;
gfx define field matrixC composite x_uv x_u x_v;

# caclulate determinants
det("matrixA","detA");
det("matrixB","detB");
det("matrixC","detC");

# define matrix based fields used in curvature equations
# (note you can now create a constant field where needed using square 
# brackets around a value, rather than having to specify it separately 
# using "gfx define field constant")
gfx define field detA_detB multiply fields detA detB;
gfx define field detC_squared power fields detC [2];

# calculate magnitudes used in curvature equations
gfx define field mag_x_u magnitude field x_u;
gfx define field mag_x_u_squared power fields mag_x_u [2];
gfx define field mag_x_v magnitude field x_v;
gfx define field mag_x_v_squared power fields mag_x_v [2];
gfx define field mag_product multiply fields mag_x_v_squared mag_x_u_squared;

# caclulate dot product based fields used in curvature equations
gfx define field dot_x_u_x_v dot_product field x_u x_v;
gfx define field dot_squared power field dot_x_u_x_v [2];

# assemble field for gaussian curvature
gfx define field gauss_numerator add fields detA_detB detC_squared scale_factors 1 -1;
gfx define field gauss_diff add fields mag_product dot_squared scale_factors 1 -1;
gfx define field gauss_denominator power field gauss_diff [2];
gfx define field gaussian_curvature divide fields gauss_numerator gauss_denominator;
gfx define field mag_gaussian_curvature magnitude field gaussian_curvature;

# assemble field for mean curvature
gfx define field mean_1 multiply fields detA mag_x_v_squared;
gfx define field mean_2 multiply fields detC dot_x_u_x_v;
gfx define field mean_lhs add fields mean_1 mean_2 scale_factors 1 -2;
gfx define field mean_rhs multiply fields detB mag_x_u_squared;
gfx define field mean_numerator add fields mean_lhs mean_rhs scale_factors 1 1;
gfx define field mean_denom power fields gauss_diff [1.5];
gfx define field mean_denominator multiply fields mean_denom [2];
gfx define field mean_curvature divide fields mean_numerator mean_denominator;
gfx define field mag_mean_curvature magnitude field mean_curvature;

# Note that the formulae for normal and gaussian curvature produces a signed
# result depending on the orientation of the surface.
# The mag_normal_curvature and mag_gaussian_curvature fields
# can be used if you want something which is always positive.

# To view the shading of the surface for mean and gaussian curvature
# either load up the scene editor and select which data is used to 
# shade the surface or use one of these commands:

#gfx modify g_element display_elements surfaces data kxi1;
#gfx modify g_element display_elements surfaces data kxi2;
#gfx modify g_element display_elements surfaces data mean_curvature_test;
#gfx modify g_element display_elements surfaces data gaussian_curvature_test;
#gfx modify g_element display_elements surfaces data mag_mean_curvature;
#gfx modify g_element display_elements surfaces data mag_gaussian_curvature;

# After using these commands you may wish to adjust your spectrum 
# using the spectrum editor or the command
#gfx modify spectrum curvature autorange;


# Similarly the values displayed at element points can be changed 
# by using the scene editor and selecting which data is used to 
# label the element points or by using of these commands:

#gfx modify g_element display_elements element_points as points glyph point label kxi1 discretization "5*2*3";
#gfx modify g_element display_elements element_points as points glyph point label kxi2 discretization "5*2*3";
#gfx modify g_element display_elements element_points as points glyph point label mean_curvature_test discretization "5*2*3";
#gfx modify g_element display_elements element_points as points glyph point label gaussian_curvature_test discretization "5*2*3";
#gfx modify g_element display_elements element_points as points glyph point label mag_mean_curvature discretization "5*2*3";
#gfx modify g_element display_elements element_points as points glyph point label mag_gaussian_curvature discretization "5*2*3";

# It is quite educational to alter the surface, and see how the curvature varies as the surface is altered.
# Use the scene editor to display nodes, using "cube_solid" as the glyph.
# On the 3D/Graphics window, use the "Node tool". Enable "edit" for this tool.
# Now drag the nodes around, and see how the curvature updates interactively.
# If you are even more adventurous, you can edit the nodal derivatives in the node file
# and see the effect. It is also possible to interactively adjust the nodal derivatives
# by defining fields for them, and then adding the appropriate glyph.
# See example_at for a demonstration.


Files used by this example are:

Name                   Modified     Size

curvature.com 17-Mar-2014 18k COPYRIGHT 17-Mar-2014 504 simple_element.exelem 17-Mar-2014 4.2k simple_element.exnode 17-Mar-2014 3.2k

Download the entire example:

Name                         Modified     Size

examples_a_curvature.tar.gz 09-Mar-2016 123k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmgui-wxFailureSun Mar 6 00:13:11 20162
last breakTue Feb 24 03:13:00 20151
cmgui-wx-debugFailureSun Mar 6 00:14:02 20162
last breakTue Feb 24 03:10:00 20152
cmgui-wx-debug-memorycheckFailureSun Mar 6 00:14:04 20162
last breakTue Feb 24 03:12:00 20153
cmgui-wx-debug-valgrindFailureSun Mar 6 01:11:23 201645
last breakSun Mar 6 01:10:00 201645
x86_64-linux
cmgui-wxFailureSun Mar 6 00:01:29 20160
last breakSun Mar 6 00:01:00 20160
cmgui-wx-debugFailureSun Mar 6 00:01:29 20160
last breakSun Mar 6 00:01:00 20160
cmgui-wx-debug-memorycheckFailureSun Mar 6 00:01:29 20160
last breakSun Mar 6 00:01:00 20160
cmgui-wx-debug-valgrindFailureSun Mar 6 00:02:41 20169
last breakSun Mar 6 00:02:00 20169
cmgui-wx-gcc-cad-debug-valgrindSuccessThu Jan 7 00:02:27 20167

Testing status by file:


Html last generated: Wed Mar 9 16:01:34 2016

Input last modified: Wed Mar 9 15:49:32 2016


CMISS Help / Examples / a / curvature