Example a_backup/mesh_integration: Calculate volume and surface area of a mesh.

This example reads a model consisting of all the standard 3-D element shapes. It defines Gauss points for a certain order Gaussian Quadrature Scheme and sums their weights multiplied by dV or dA (differential volume or area at the Gauss point) to calculate the total volume and surface area of the mesh.

Screenshot of example a_backup/mesh_integration


The comfile run by this example is as follows:

# ***** BEGIN LICENSE BLOCK *****
# Version: MPL 1.1/GPL 2.0/LGPL 2.1
#
# The contents of this file are subject to the Mozilla Public License Version
# 1.1 (the "License"); you may not use this file except in compliance with
# the License. You may obtain a copy of the License at
# http://www.mozilla.org/MPL/
#
# Software distributed under the License is distributed on an "AS IS" basis,
# WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
# for the specific language governing rights and limitations under the
# License.
#
# The Original Code is CMISS Mesh Integration example.
#
# The Initial Developer of the Original Code is
# Auckland Uniservices Ltd, Auckland, New Zealand.
# Portions created by the Initial Developer are Copyright (C) 2011
# the Initial Developer. All Rights Reserved.
#
# Alternatively, the contents of this file may be used under the terms of
# either the GNU General Public License Version 2 or later (the "GPL"), or
# the GNU Lesser General Public License Version 2.1 or later (the "LGPL"),
# in which case the provisions of the GPL or the LGPL are applicable instead
# of those above. If you wish to allow use of your version of this file only
# under the terms of either the GPL or the LGPL, and not to allow others to
# use your version of this file under the terms of the MPL, indicate your
# decision by deleting the provisions above and replace them with the notice
# and other provisions required by the GPL or the LGPL. If you do not delete
# the provisions above, a recipient may use your version of this file under
# the terms of any one of the MPL, the GPL or the LGPL.
#
# ***** END LICENSE BLOCK *****

# This example calculates volume and surface area of a mesh using Gauss quadrature.
# It requires cmgui version 2.8 (October 2011) or later.

# Read mesh containing all possible 3-D element shapes
gfx read nodes $example/allshapes.exnode;
gfx read elements $example/allshapes.exelem;

# Create an element group containing just the exterior faces, and clear graphics for it
gfx create egroup exterior region allshapes;
gfx modify egroup allshapes/exterior add faces 2,4..5,8..10,12..14,16..24;
gfx modify g_element allshapes/exterior general clear;

# create order 2 Gauss points on both the volume and exterior surface to integrate with
gfx define field allshapes/gauss_location finite_element number_of_components 1 field element_xi;
gfx define field allshapes/gauss_weight finite_element number_of_components 1 field real;
gfx create dgroup gausspoints region allshapes;
gfx create gauss_points region allshapes mesh cmiss_mesh_3d order 2 gauss_location_field gauss_location gauss_point_nodeset gausspoints.cmiss_data gauss_weight_field gauss_weight first_identifier 1;
gfx create dgroup surface_gausspoints region allshapes;
gfx create gauss_points region allshapes mesh exterior.cmiss_mesh_2d order 2 gauss_location_field gauss_location gauss_point_nodeset surface_gausspoints.cmiss_data gauss_weight_field gauss_weight first_identifier 10001;

# get derivatives of coordinate field w.r.t. element chart 'xi'
# (Be aware that the 'xi' is for the element you evaluate it at, so differs for 3-D and 2-D)
gfx define field allshapes/dX_dxi1 basis_derivative fe_field coordinates order 1 xi_indices 1;
gfx define field allshapes/dX_dxi2 basis_derivative fe_field coordinates order 1 xi_indices 2;
gfx define field allshapes/dX_dxi3 basis_derivative fe_field coordinates order 1 xi_indices 3;
# pack these into a 3x3 matrix, and get differential volume dV from its determinant, in 3-D elements
gfx define field allshapes/dX_dxi composite dX_dxi1 dX_dxi2 dX_dxi3;
gfx define field allshapes/dV determinant field dX_dxi;
# magnitude of dX_dxi1 (x) dX_dxi2 gives dA when evaluated on 2-D elements
gfx define field allshapes/dX_dxi_cross cross_product fields dX_dxi1 dX_dxi2;
gfx define field allshapes/dA magnitude field dX_dxi_cross;

# obtain fields at gauss points:
gfx define field allshapes/gauss_coordinates embedded element_xi gauss_location field coordinates;
gfx define field allshapes/gauss_dV embedded element_xi gauss_location field dV;
gfx define field allshapes/gauss_dA embedded element_xi gauss_location field dA;

# sum weighted dV over gauss points
gfx define field allshapes/weighted_gauss_dV multiply fields gauss_weight gauss_dV;
gfx define field allshapes/volume nodeset_sum field weighted_gauss_dV nodeset gausspoints.cmiss_data;

# sum weighted dA over surface gauss points
gfx define field allshapes/weighted_gauss_dA multiply fields gauss_weight gauss_dA;
gfx define field allshapes/area nodeset_sum field weighted_gauss_dA nodeset surface_gausspoints.cmiss_data;

# visualise mesh and both sets of Gauss points
gfx modify g_element allshapes/gausspoints general clear;
gfx modify g_element allshapes/gausspoints data_points coordinate gauss_coordinates LOCAL glyph cross general size "0.05*0.05*0.05" centre 0,0,0 font default select_on material silver selected_material default_selected;
gfx modify g_element allshapes/surface_gausspoints general clear;
gfx modify g_element allshapes/surface_gausspoints data_points coordinate gauss_coordinates LOCAL glyph cross general size "0.05*0.05*0.05" centre 0,0,0 font default select_on material green selected_material default_selected;
gfx modify g_element allshapes general clear;
gfx modify g_element allshapes lines coordinate coordinates tessellation default LOCAL native_discretization NONE select_on material default selected_material default_selected;
gfx modify g_element allshapes node_points coordinate coordinates LOCAL glyph sphere general size "0.07*0.07*0.07" centre 0,0,0 font default select_on material gold selected_material default_selected;

# show mesh volume in silver in the window corner 
gfx modify g_element allshapes point coordinate coordinates NORMALISED_WINDOW_FIT_LEFT glyph none general size "1*1*1" centre 0.95,-0.9,0 font default label volume select_on material silver selected_material default_selected;
# show mesh surface area in silver in the window corner
gfx modify g_element allshapes point coordinate coordinates NORMALISED_WINDOW_FIT_LEFT glyph none general size "1*1*1" centre 0.95,-0.8,0 font default label area select_on material green selected_material default_selected;

if (!$TESTING) {
	gfx create window 1;
}

# Clear all the gauss points and recreate for order 3 quadrature:
gfx destroy data all group allshapes;
gfx create gauss_points region allshapes mesh cmiss_mesh_3d order 3 gauss_location_field gauss_location gauss_point_nodeset gausspoints.cmiss_data gauss_weight_field gauss_weight first_identifier 1;
gfx create gauss_points region allshapes mesh exterior.cmiss_mesh_2d order 3 gauss_location_field gauss_location gauss_point_nodeset surface_gausspoints.cmiss_data gauss_weight_field gauss_weight first_identifier 10001;

# It is instructive to inspect the fields that are summed here.
# Show field gauss_weight, gauss_location, gauss_dA and gauss_dV as labels on the data_points
# Also try editing the nodes using the node tool in the graphics window
# to see the volume and area change

if ($TESTING) {
	# Test all quadrature orders produce the same volume and area
	my $scaling;
	my $prefix = "";
	for ($scale = 1; $scale <= 2; $scale++) {
		my $order;
		for ($order = 1; $order <= 4; $order++) {
			gfx destroy data all group allshapes;
			gfx create gauss_points region allshapes mesh cmiss_mesh_3d order $order gauss_location_field gauss_location gauss_point_nodeset gausspoints.cmiss_data gauss_weight_field gauss_weight first_identifier 1;
			gfx create gauss_points region allshapes mesh exterior.cmiss_mesh_2d order $order gauss_location_field gauss_location gauss_point_nodeset surface_gausspoints.cmiss_data gauss_weight_field gauss_weight first_identifier 10001;

			# for testing answers use trick of evaluating into constant field, wastefully for all nodes

			my $volume_name = $prefix . "volume_order$order";
			gfx define field "allshapes/$volume_name" constant 1;
			gfx evaluate ngroup allshapes source volume destination $volume_name;
			gfx list field "allshapes/$volume_name" commands;

			my $area_name = $prefix . "area_order$order";
			gfx define field "allshapes/$area_name" constant 1;
			gfx evaluate ngroup allshapes source area destination $area_name;
			gfx list field "allshapes/$area_name" commands;
		}
		if ($scale == 1) {
			$prefix = "scaled_";
			gfx define field allshapes/scaled_coordinates scale field coordinates scale_factors 1.5 0.75 2.0;
			gfx evaluate ngroup allshapes source scaled_coordinates destination coordinates;
		}
	}
	# restore original nodal coordinates
	gfx read nodes $example/allshapes.exnode;
}

Files used by this example are:

Name                  Modified     Size

mesh_integration.com 20-Apr-2012 8.5k COPYRIGHT 20-Apr-2012 504 allshapes.exelem 20-Apr-2012 10k allshapes.exnode 20-Apr-2012 1.1k

Download the entire example:

Name                                       Modified     Size

examples_a_backup_mesh_integration.tar.gz 12-Aug-2014 24k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmgui-wxSuccessSun Mar 6 00:06:34 20160
cmgui-wx-debugSuccessSun Mar 6 00:06:35 20161
cmgui-wx-debug-memorycheckSuccessSun Mar 6 00:06:35 20161
cmgui-wx-debug-valgrindSuccessSun Mar 6 00:39:18 201624
x86_64-linux
cmgui-wxSuccessThu Jan 7 00:01:29 20160
cmgui-wx-debugSuccessThu Jan 7 00:01:29 20160
cmgui-wx-debug-memorycheckSuccessThu Jan 7 00:01:29 20160
cmgui-wx-debug-valgrindFailureSun Mar 6 00:04:24 201621
last breakSun Mar 6 00:04:00 201621
last successWed Jun 3 00:18:00 201520
cmgui-wx-gcc-cad-debug-valgrindSuccessThu Jan 7 00:06:35 201629

Testing status by file:


Html last generated: Sun Mar 6 05:50:59 2016

Input last modified: Fri Apr 20 16:00:17 2012


CMISS Help / Examples / a_backup / mesh_integration