[cig-commits] r5160 - in long/3D/Gale/trunk/src/Underworld: .
Utils/src
walter at geodynamics.org
walter at geodynamics.org
Tue Oct 31 13:33:23 PST 2006
Author: walter
Date: 2006-10-31 13:33:22 -0800 (Tue, 31 Oct 2006)
New Revision: 5160
Added:
long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.c
long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.h
long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.meta
Modified:
long/3D/Gale/trunk/src/Underworld/
long/3D/Gale/trunk/src/Underworld/Utils/src/Init.c
long/3D/Gale/trunk/src/Underworld/Utils/src/Utils.h
long/3D/Gale/trunk/src/Underworld/Utils/src/types.h
Log:
r698 at earth: boo | 2006-10-31 13:32:08 -0800
r688 at earth (orig r365): MirkoVelic | 2006-10-30 07:58:26 -0800
The Voronoi integration scheme.
This needs to be more thoroughly tested.
The 3D mode seems to work fine for one model
at least.
Property changes on: long/3D/Gale/trunk/src/Underworld
___________________________________________________________________
Name: svk:merge
- 9570c393-cf10-0410-b476-9a651db1e55a:/cig:697
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:364
+ 9570c393-cf10-0410-b476-9a651db1e55a:/cig:698
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:365
Added: long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.c 2006-10-31 21:33:19 UTC (rev 5159)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.c 2006-10-31 21:33:22 UTC (rev 5160)
@@ -0,0 +1,926 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+** Copyright (c) 2006, Monash Cluster Computing
+** All rights reserved.
+** Redistribution and use in source and binary forms, with or without modification,
+** are permitted provided that the following conditions are met:
+**
+** * Redistributions of source code must retain the above copyright notice,
+** this list of conditions and the following disclaimer.
+** * Redistributions in binary form must reproduce the above copyright
+** notice, this list of conditions and the following disclaimer in the
+** documentation and/or other materials provided with the distribution.
+** * Neither the name of the Monash University nor the names of its contributors
+** may be used to endorse or promote products derived from this software
+** without specific prior written permission.
+**
+** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
+** THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+** PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
+** BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+** OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+**
+**
+** Contact:
+*% Louis Moresi - Louis.Moresi at sci.monash.edu.au
+*%
+** Author:
+** Mirko Velic - Mirko.Velic at sci.monash.edu.au
+**
+** Assumptions:
+** I am assuming that the xi's (local coords) on the IntegrationPoint particles
+** are precalculated somewhere and get reset based on material PIC positions each time step.
+**
+** Notes:
+** The DVCWeights class should really be a class the next level up here.
+** We should be able to swap out the WeightsCalculator_CalculateAll function instead of just setting
+** a pointer inside that function
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+
+/****************************************************************************************************************
+
+ The algorithm here-in computes a discrete voronoi diagram per FEM cell given a set of local
+ particle positions, in 3D and 2D. The volumes of the Voronoi regions are used as integration weights for
+ the integration point swarm and the integration points are the centroids of the same volumes.
+
+ For a description of this algorithm, see the article by Velic et.al.
+ "A Fast Robust Algorithm for computing Discrete Voronoi Diagrams in N-dimensions"
+
+ Warning:
+
+ Be very careful of making changes here. It may have undesirable consequences in regard to the
+ speed of this implementation. Speed is very important here. You may be tempted to merge the
+ 2D and 3D functions. Don't do it unless you can do it without using control statements.
+ This implementation uses von-Neumann neighbourhoods for boundary chain growth. Do not be tempted
+ to implement "diagonal" neighbourhood growth cycles. For this is an abomination unto me.
+*****************************************************************************************************************/
+
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StG_FEM/StG_FEM.h>
+
+#include "types.h"
+#include <PICellerator/PICellerator.h>
+#include "DVCWeights.h"
+
+#include <assert.h>
+#include <string.h>
+#include <math.h>
+
+
+/* Textual name of this class */
+const Type DVCWeights_Type = "DVCWeights";
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Constructors
+*/
+DVCWeights* DVCWeights_New( Name name, Dimension_Index dim ) {
+ DVCWeights* self = (DVCWeights*) _DVCWeights_DefaultNew( name );
+
+ DVCWeights_InitAll( self, dim );
+ return self;
+}
+
+DVCWeights* _DVCWeights_New(
+ SizeT _sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ WeightsCalculator_CalculateFunction* _calculate,
+ Name name )
+{
+ DVCWeights* self;
+
+ /* Allocate memory */
+ assert( _sizeOfSelf >= sizeof(DVCWeights) );
+ self = (DVCWeights*)_WeightsCalculator_New(
+ _sizeOfSelf,
+ type,
+ _delete,
+ _print,
+ _copy,
+ _defaultConstructor,
+ _construct,
+ _build,
+ _initialise,
+ _execute,
+ _destroy,
+ _calculate,
+ name );
+
+
+ /* General info */
+
+ /* Virtual Info */
+ return self;
+}
+
+void _DVCWeights_Init( void* dvcWeights, int *res ) {
+ DVCWeights* self = (DVCWeights*)dvcWeights;
+ self->isConstructed = True;
+
+ self->resX = res[I_AXIS];
+ self->resY = res[J_AXIS];
+ self->resZ = res[K_AXIS];
+
+}
+
+void DVCWeights_InitAll( void* dvcWeights, Dimension_Index dim ) {
+ DVCWeights* self = (DVCWeights*)dvcWeights;
+ WeightsCalculator_InitAll( self, dim );
+}
+
+/*------------------------------------------------------------------------------------------------------------------------
+** Virtual functions
+*/
+
+void _DVCWeights_Delete( void* dvcWeights ) {
+ DVCWeights* self = (DVCWeights*)dvcWeights;
+ /* Delete parent */
+ _WeightsCalculator_Delete( self );
+}
+
+
+void _DVCWeights_Print( void* dvcWeights, Stream* stream ) {
+ DVCWeights* self = (DVCWeights*)dvcWeights;
+ /* Print parent */
+ _WeightsCalculator_Print( self, stream );
+}
+
+
+
+void* _DVCWeights_Copy( void* dvcWeights, void* dest, Bool deep, Name nameExt, PtrMap* ptrMap ) {
+ DVCWeights* self = (DVCWeights*)dvcWeights;
+ DVCWeights* newDVCWeights;
+
+ newDVCWeights = (DVCWeights*)_WeightsCalculator_Copy( self, dest, deep, nameExt, ptrMap );
+ return (void*)newDVCWeights;
+}
+
+void* _DVCWeights_DefaultNew( Name name ) {
+ return (void*) _DVCWeights_New(
+ sizeof(DVCWeights),
+ DVCWeights_Type,
+ _DVCWeights_Delete,
+ _DVCWeights_Print,
+ _DVCWeights_Copy,
+ _DVCWeights_DefaultNew,
+ _DVCWeights_Construct,
+ _DVCWeights_Build,
+ _DVCWeights_Initialise,
+ _DVCWeights_Execute,
+ _DVCWeights_Destroy,
+ _DVCWeights_Calculate,
+ name );
+}
+
+
+void _DVCWeights_Construct( void* dvcWeights, Stg_ComponentFactory* cf ) {
+ DVCWeights* self = (DVCWeights*) dvcWeights;
+
+ int defaultResolution;
+ int resolution[3];
+
+
+ defaultResolution = Stg_ComponentFactory_GetUnsignedInt( cf, self->name, "resolution", 10 );
+ resolution[ I_AXIS ] = Stg_ComponentFactory_GetUnsignedInt( cf, self->name, "resolutionX", defaultResolution );
+ resolution[ J_AXIS ] = Stg_ComponentFactory_GetUnsignedInt( cf, self->name, "resolutionY", defaultResolution );
+ resolution[ K_AXIS ] = Stg_ComponentFactory_GetUnsignedInt( cf, self->name, "resolutionZ", defaultResolution );
+
+ _WeightsCalculator_Construct( self, cf );
+
+ _DVCWeights_Init( self, resolution );
+}
+
+void _DVCWeights_Build( void* dvcWeights, void* data ) {
+ DVCWeights* self = (DVCWeights*)dvcWeights;
+ _WeightsCalculator_Build( self, data );
+}
+void _DVCWeights_Initialise( void* dvcWeights, void* data ) {
+ DVCWeights* self = (DVCWeights*)dvcWeights;
+ _WeightsCalculator_Initialise( self, data );
+}
+void _DVCWeights_Execute( void* dvcWeights, void* data ) {
+ DVCWeights* self = (DVCWeights*)dvcWeights;
+ _WeightsCalculator_Execute( self, data );
+}
+void _DVCWeights_Destroy( void* dvcWeights, void* data ) {
+ DVCWeights* self = (DVCWeights*)dvcWeights;
+ _WeightsCalculator_Destroy( self, data );
+}
+
+/*-------------------------------------------------------------------------------------------------------------------------
+** Private Functions
+*/
+
+void get_centroids( struct cell *cells,struct particle *pList,
+ int n, int m, int l,int nump,double vol){
+ int i;
+ int *count;
+
+ count=(int *)malloc(sizeof(int)*nump);
+
+ for(i=0;i<nump;i++){
+ count[i] = 0;
+ pList[ i ].cx = 0.0;
+ pList[ i ].cy = 0.0;
+ pList[ i ].cz = 0.0;
+ }
+ for(i=0;i<n*m*l;i++){
+ pList[ cells[i].p ].cx += cells[i].x;
+ pList[ cells[i].p ].cy += cells[i].y;
+ pList[ cells[i].p ].cz += cells[i].z;
+ count[ cells[i].p ]++;//for total volume of a cell
+ }
+ for(i=0;i<nump;i++){
+ pList[ i ].w = count[i]*vol;
+ if(count[i] != 0){
+ pList[ i ].cx /= count[i];
+ pList[ i ].cy /= count[i];
+ pList[ i ].cz /= count[i];
+ }
+ }
+ free(count);
+}
+void get_centroids2d( struct cell2d *cells,struct particle2d *pList,
+ int n, int m, int nump,double vol){
+ int i;
+ int *count;
+
+ count=(int *)malloc(sizeof(int)*nump);
+
+ for(i=0;i<nump;i++){
+ count[i] = 0;
+ pList[ i ].cx = 0.0;
+ pList[ i ].cy = 0.0;
+ }
+ for(i=0;i<n*m;i++){
+ pList[ cells[i].p ].cx += cells[i].x;
+ pList[ cells[i].p ].cy += cells[i].y;
+ count[ cells[i].p ]++;//for total volume of a cell
+ }
+ for(i=0;i<nump;i++){
+ pList[ i ].w = count[i]*vol;
+ if(count[i] != 0){
+ pList[ i ].cx /= count[i];
+ pList[ i ].cy /= count[i];
+ }
+ }
+ free(count);
+}
+
+void claim_cells(struct chain **bbchain,struct cell **ccells,struct particle **pList,int p_i){
+ int i,count;
+ int cell_num0;
+ double x0,y0,x1,y1,x2,y2,z0,z1,z2,dist1,dist2;
+ struct chain *bchain = &(*bbchain)[p_i];
+ struct cell *cells = *ccells;
+ int *temp;
+
+ count = 0;
+ bchain->numclaimed = 0;
+
+
+ for(i=0;i<bchain->sizeofboundary;i++){
+ cell_num0 = bchain->new_bound_cells[i];// cell number we are trying to claim
+ if(cells[cell_num0].p == -1){//if cell unowned then claim cell
+ /* This is the bit needed for mallocing */
+ /* do a test here to see if we need to realloc bchain->new_claimed_cells and bchain->new_bound_cells */
+ if( count > bchain->new_claimed_cells_malloced - 1 ){
+ temp = (int *)realloc( bchain->new_claimed_cells, (bchain->new_claimed_cells_malloced + INC)*sizeof(int) );
+ bchain->new_claimed_cells = temp;
+ bchain->new_claimed_cells_malloced += INC;
+ temp = (int *)realloc( bchain->new_bound_cells, (bchain->new_bound_cells_malloced + INC)*sizeof(int) );
+ bchain->new_bound_cells = temp;
+ bchain->new_bound_cells_malloced += INC;
+ }
+ /* end of bit needed for mallocing */
+ bchain->new_claimed_cells[count] = cell_num0;
+ bchain->numclaimed++;
+ count++;
+ cells[cell_num0].p = p_i;// this cell is now owned by particle p_i
+ }
+ else{
+ if(cells[cell_num0].p != p_i){
+ //we need a contest between particles for the cell.
+ x2 = (*pList)[p_i].x;
+ y2 = (*pList)[p_i].y;
+ z2 = (*pList)[p_i].z;
+ x1 = (*pList)[cells[cell_num0].p].x;
+ y1 = (*pList)[cells[cell_num0].p].y;
+ z1 = (*pList)[cells[cell_num0].p].z;
+ x0 = cells[cell_num0].x;
+ y0 = cells[cell_num0].y;
+ z0 = cells[cell_num0].z;
+
+ dist1 = distance(x0,y0,z0,x1,y1,z1);
+ dist2 = distance(x0,y0,z0,x2,y2,z2);
+ if(dist1 > dist2){
+ bchain->new_claimed_cells[count] = cell_num0;
+ bchain->numclaimed++;
+ count++;
+ cells[cell_num0].p = p_i;// this cell is now owned by particle p_i
+ }
+ }//if
+ }//else
+ }
+ bchain->new_claimed_cells[count] = -1;// end of list
+}
+
+void claim_cells2d(struct chain **bbchain,struct cell2d **ccells,struct particle2d **pList,int p_i){
+ int i,count;
+ int cell_num0;
+ double x0,y0,x1,y1,x2,y2,dist1,dist2;
+ struct chain *bchain = &(*bbchain)[p_i];
+ struct cell2d *cells = *ccells;
+ int *temp;
+
+ count = 0;
+ bchain->numclaimed = 0;
+
+
+ for(i=0;i<bchain->sizeofboundary;i++){
+ cell_num0 = bchain->new_bound_cells[i];// cell number we are trying to claim
+ if(cells[cell_num0].p == -1){//if cell unowned then claim cell
+ /* This is the bit needed for mallocing */
+ /* do a test here to see if we need to realloc bchain->new_claimed_cells and bchain->new_bound_cells */
+ if( count > bchain->new_claimed_cells_malloced - 1 ){
+ temp = (int *)realloc( bchain->new_claimed_cells, (bchain->new_claimed_cells_malloced + INC)*sizeof(int) );
+ bchain->new_claimed_cells = temp;
+ bchain->new_claimed_cells_malloced += INC;
+ temp = (int *)realloc( bchain->new_bound_cells, (bchain->new_bound_cells_malloced + INC)*sizeof(int) );
+ bchain->new_bound_cells = temp;
+ bchain->new_bound_cells_malloced += INC;
+ }
+ /* end of bit needed for mallocing */
+ bchain->new_claimed_cells[count] = cell_num0;
+ bchain->numclaimed++;
+ count++;
+ cells[cell_num0].p = p_i;// this cell is now owned by particle p_i
+ }
+ else{
+ if(cells[cell_num0].p != p_i){
+ //we need a contest between particles for the cell.
+ x2 = (*pList)[p_i].x;
+ y2 = (*pList)[p_i].y;
+ x1 = (*pList)[cells[cell_num0].p].x;
+ y1 = (*pList)[cells[cell_num0].p].y;
+ x0 = cells[cell_num0].x;
+ y0 = cells[cell_num0].y;
+
+ dist1 = distance2d(x0,y0,x1,y1);
+ dist2 = distance2d(x0,y0,x2,y2);
+ if(dist1 > dist2){
+ bchain->new_claimed_cells[count] = cell_num0;
+ bchain->numclaimed++;
+ count++;
+ cells[cell_num0].p = p_i;// this cell is now owned by particle p_i
+ }
+ }//if
+ }//else
+ }
+ bchain->new_claimed_cells[count] = -1;// end of list
+}
+
+
+void reset_grid(struct cell **cells, int n){
+ int i;
+
+ for(i=0;i<n;i++){
+ (*cells)[i].p = -1;
+ (*cells)[i].done = 0;
+ }
+}
+void reset_grid2d(struct cell2d **cells, int n){
+ int i;
+
+ for(i=0;i<n;i++){
+ (*cells)[i].p = -1;
+ (*cells)[i].done = 0;
+ }
+}
+void update_bchain(struct chain **bbchain,struct cell **ccells,int p_i){
+ int i,k,count;
+ int cell_num0,cell_num[6],cell_num1;
+ struct chain *bchain = &(*bbchain)[p_i];
+ struct cell *cells = *ccells;
+ int *temp;
+
+ count = 0;
+ bchain->sizeofboundary = 0;
+ for(i=0;i<bchain->numclaimed;i++){
+ cell_num0 =bchain->new_claimed_cells[i];
+
+ cell_num[0] = cells[cell_num0].S;
+ cell_num[1] = cells[cell_num0].N;
+ cell_num[2] = cells[cell_num0].E;
+ cell_num[3] = cells[cell_num0].W;
+ cell_num[4] = cells[cell_num0].U;
+ cell_num[5] = cells[cell_num0].D;
+
+ for(k=0;k<6;k++){
+ cell_num1 = cell_num[k];
+ // if cell does not already belong to the particle and hasn't been
+ // marked as being done then add it to new boundary array and mark it
+ // as done
+ if(cell_num1 != -2){
+ if(cells[cell_num1].p != p_i && cells[cell_num1].done != 1){
+ /* This is the bit needed for mallocing */
+ /* do a test here to see if we need to realloc bchain->new_claimed_cells and bchain->new_bound_cells */
+ if( count > bchain->new_bound_cells_malloced - 1 ){
+ temp = (int *)realloc( bchain->new_claimed_cells, (bchain->new_claimed_cells_malloced + INC)*sizeof(int) );
+ bchain->new_claimed_cells = temp;
+ bchain->new_claimed_cells_malloced += INC;
+ temp = (int *)realloc( bchain->new_bound_cells, (bchain->new_bound_cells_malloced + INC)*sizeof(int) );
+ bchain->new_bound_cells = temp;
+ bchain->new_bound_cells_malloced += INC;
+ }
+ /* end of bit needed for mallocing */
+ bchain->new_bound_cells[count] = cell_num1;
+ bchain->sizeofboundary++;
+ count++;
+ cells[cell_num1].done = 1;
+ }//if
+ }//if cell_num1
+ }//for k
+ }//for
+ // reset the done flags back to zero for next time
+ for(i=0;i<count;i++){
+ cells[ bchain->new_bound_cells[i] ].done = 0;
+ }
+}
+
+void update_bchain2d(struct chain **bbchain,struct cell2d **ccells,int p_i){
+ int i,k,count;
+ int cell_num0,cell_num[4],cell_num1;
+ struct chain *bchain = &(*bbchain)[p_i];
+ struct cell2d *cells = *ccells;
+ int *temp;
+
+ count = 0;
+ bchain->sizeofboundary = 0;
+ for(i=0;i<bchain->numclaimed;i++){
+ cell_num0 =bchain->new_claimed_cells[i];
+
+ cell_num[0] = cells[cell_num0].S;
+ cell_num[1] = cells[cell_num0].N;
+ cell_num[2] = cells[cell_num0].E;
+ cell_num[3] = cells[cell_num0].W;
+
+ for(k=0;k<4;k++){
+ cell_num1 = cell_num[k];
+ // if cell does not already belong to the particle and hasn't been
+ // marked as being done then add it to new boundary array and mark it
+ // as done
+ if(cell_num1 != -2){
+ if(cells[cell_num1].p != p_i && cells[cell_num1].done != 1){
+ /* This is the bit needed for mallocing */
+ /* do a test here to see if we need to realloc bchain->new_claimed_cells and bchain->new_bound_cells */
+ if( count > bchain->new_bound_cells_malloced - 1 ){
+ temp = (int *)realloc( bchain->new_claimed_cells, (bchain->new_claimed_cells_malloced + INC)*sizeof(int) );
+ bchain->new_claimed_cells = temp;
+ bchain->new_claimed_cells_malloced += INC;
+ temp = (int *)realloc( bchain->new_bound_cells, (bchain->new_bound_cells_malloced + INC)*sizeof(int) );
+ bchain->new_bound_cells = temp;
+ bchain->new_bound_cells_malloced += INC;
+ }
+ /* end of bit needed for mallocing */
+ bchain->new_bound_cells[count] = cell_num1;
+ bchain->sizeofboundary++;
+ count++;
+ cells[cell_num1].done = 1;
+ }//if
+ }//if cell_num1
+ }//for k
+ }//for
+ // reset the done flags back to zero for next time
+ for(i=0;i<count;i++){
+ cells[ bchain->new_bound_cells[i] ].done = 0;
+ }
+}
+
+
+/***********************************************************
+ n m l are z x y grid resolutions respectively.
+************************************************************/
+void construct_grid(struct cell **cell_list, int n, int m, int l,
+ double x0,double y0,double z0,double x1,double y1,double z1){
+ struct cell *cells;
+ int i,j,k;
+ double dx,dy,dz,Dx,Dy,Dz,X,Y,Z;
+
+ cells = malloc(n*m*l*sizeof(struct cell));
+ for(i=0;i<l*m*n;i++){
+ cells[i].S = -2;
+ cells[i].N = -2;
+ cells[i].E = -2;
+ cells[i].W = -2;
+ cells[i].U = -2;
+ cells[i].D = -2;
+ cells[i].p = -1;
+ cells[i].done = 0;
+ }
+ for(k=0;k<n;k++){
+ for(i=0;i<l*(m-1);i++){
+ cells[i+k*l*m].N=i+k*l*m+l;
+ cells[i+l+k*l*m].S = i+k*l*m;
+ }
+ }
+ for(k=0;k<n-1;k++){
+ for(i=0;i<l*m;i++){
+ cells[i+k*l*m].U=i+k*l*m+l*m;
+ cells[i+k*l*m+l*m].D = i+k*l*m;
+ }
+ }
+
+ Dx = x1-x0; Dy = y1-y0; Dz = z1-z0;
+ dx = Dx/l; dy = Dy/m; dz = Dz/n;
+ Z = z0 - dz/2.0;
+ for(k=0;k<n;k++){
+ Z = Z + dz;
+ Y = y0 - dy/2.0;
+ for(j=0;j<m;j++){
+ Y = Y + dy;
+ X = x0 - dx/2.0;
+ for(i=0;i<l;i++){
+ X=X+dx;
+ cells[i+l*j+k*l*m].x = X;
+ cells[i+l*j+k*l*m].y = Y;
+ cells[i+l*j+k*l*m].z = Z;
+ if(i!= l-1){
+ cells[i+l*j+k*l*m].E = i+l*j+1+k*l*m;
+ cells[i+l*j+1+k*l*m].W = i+l*j+k*l*m;
+ }
+ }//x
+ }//y
+ }//z
+ *cell_list = cells;
+}
+
+/***********************************************************
+ m l are x y grid resolutions respectively.
+************************************************************/
+void construct_grid2d(struct cell2d **cell_list, int m, int l,
+ double x0,double y0,double x1,double y1){
+ struct cell2d *cells;
+ int i,j;
+ double dx,dy,Dx,Dy,X,Y;
+
+ cells = malloc(m*l*sizeof(struct cell));
+ for(i=0;i<l*m;i++){
+ cells[i].S = -2;
+ cells[i].N = -2;
+ cells[i].E = -2;
+ cells[i].W = -2;
+ cells[i].p = -1;
+ cells[i].done = 0;
+ }
+ for(i=0;i<l*(m-1);i++){
+ cells[i].N=i+l;
+ cells[i+l].S = i;
+ }
+ Dx = x1-x0; Dy = y1-y0;
+ dx = Dx/l; dy = Dy/m;
+ Y = y0 - dy/2.0;
+ for(j=0;j<m;j++){
+ Y = Y + dy;
+ X = x0 - dx/2.0;
+ for(i=0;i<l;i++){
+ X=X+dx;
+ cells[i+l*j].x = X;
+ cells[i+l*j].y = Y;
+ if(i!= l-1){
+ cells[i+l*j].E = i+l*j+1;
+ cells[i+l*j+1].W = i+l*j;
+ }
+ }
+ }
+ *cell_list = cells;
+}
+
+
+double distance(double x0, double y0, double z0, double x1, double y1, double z1){
+ return (x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)+(z1-z0)*(z1-z0);
+}
+double distance2d(double x0, double y0,double x1, double y1){
+ return (x1-x0)*(x1-x0)+(y1-y0)*(y1-y0);
+}
+
+void init_structs( struct chain **bchain, struct particle **pList, int nump){
+
+ int i;
+ // init the data structures
+ if( (*bchain = (struct chain *)malloc( nump*sizeof(struct chain ) )) == 0){
+ printf("No memory for bchain\nCannot continue.\n");
+ exit(1);
+ }
+ // note that doing bchain[i]->new... doesn't work
+ for(i=0;i<nump;i++){
+ (*bchain)[i].new_claimed_cells = (int *)malloc(INC*sizeof(int));
+ (*bchain)[i].new_claimed_cells_malloced = INC;
+ (*bchain)[i].new_bound_cells = (int *)malloc(INC*sizeof(int));
+ (*bchain)[i].new_bound_cells_malloced = INC;
+ }
+ if( (*pList = (struct particle *)malloc( nump*sizeof(struct particle ) )) == 0){
+ printf("No memory for pList\nCannot continue.\n");
+ exit(1);
+ }
+
+}
+
+void init_structs2d( struct chain **bchain, struct particle2d **pList, int nump){
+
+ int i;
+ // init the data structures
+ if( (*bchain = (struct chain *)malloc( nump*sizeof(struct chain ) )) == 0){
+ printf("No memory for bchain\nCannot continue.\n");
+ exit(1);
+ }
+ //
+ for(i=0;i<nump;i++){
+ (*bchain)[i].new_claimed_cells = (int *)malloc(INC*sizeof(int));
+ (*bchain)[i].new_claimed_cells_malloced = INC;
+ (*bchain)[i].new_bound_cells = (int *)malloc(INC*sizeof(int));
+ (*bchain)[i].new_bound_cells_malloced = INC;
+ }
+ if( (*pList = (struct particle2d *)malloc( nump*sizeof(struct particle2d ) )) == 0){
+ printf("No memory for pList\nCannot continue.\n");
+ exit(1);
+ }
+
+}
+/********************************************************************
+
+ All the parameters passed into create_voronoi
+ must be initialised already.
+
+*********************************************************************/
+void create_voronoi( struct chain **bchain, struct particle **pList, struct cell **cells,
+ double dx, double dy, double dz,
+ int nump,
+ int numx, int numy, int numz,
+ double BBXMIN, double BBXMAX,
+ double BBYMIN, double BBYMAX,
+ double BBZMIN, double BBZMAX){
+ int i,j,k,l;
+ int count;
+ int claimed;
+
+ for(i=0;i<nump;i++){
+ k = ((*pList)[i].x - BBXMIN)/dx;
+ j = ((*pList)[i].y - BBYMIN)/dy;
+ l = ((*pList)[i].z - BBZMIN)/dz;
+ (*cells)[k+j*numx+l*numx*numy].p = i; //particle number i
+ (*bchain)[i].numclaimed = 1;// number of most recently claimed cells
+ (*bchain)[i].sizeofboundary = 0;
+ (*bchain)[i].totalclaimed = 1;// total of claimed cells so far.
+ (*bchain)[i].done = 0;
+ (*bchain)[i].index = k+j*numx+l*numx*numy;// ith particle is in cell # k+j*numx
+ (*bchain)[i].new_claimed_cells[0] = k+j*numx+l*numx*numy;
+ // ith particle has just claimed cell number k+j*numx+l*numx*numy
+ (*bchain)[i].new_claimed_cells[1] = -1;// denotes end of claimed_cells list
+ // when we have finished claiming cells we call this function.
+ update_bchain(bchain,cells,i);
+ }
+ count = i;// number of particles
+ claimed = 1;
+
+ while(claimed != 0){
+ claimed = 0 ;
+ for(i=0;i<count;i++){
+ claim_cells(bchain,cells,pList,i);
+ claimed += (*bchain)[i].numclaimed;
+ update_bchain(bchain,cells,i);
+ }
+ }//while
+}
+/********************************************************************
+
+ All the parameters passed into create_voronoi
+ must be initialised already.
+
+*********************************************************************/
+void create_voronoi2d( struct chain **bchain, struct particle2d **pList, struct cell2d **cells,
+ double dx, double dy,
+ int nump,
+ int numx, int numy,
+ double BBXMIN, double BBXMAX,
+ double BBYMIN, double BBYMAX){
+ int i,j,k;
+ int claimed;
+
+ for(i=0;i<nump;i++){
+ k = ((*pList)[i].x - BBXMIN)/dx;
+ j = ((*pList)[i].y - BBYMIN)/dy;
+
+ (*cells)[k+j*numx].p = i; //particle number i
+
+ (*bchain)[i].numclaimed = 1;// number of most recently claimed cells
+ (*bchain)[i].sizeofboundary = 0;
+ (*bchain)[i].totalclaimed = 1;// total of claimed cells so far.
+ (*bchain)[i].done = 0;
+ (*bchain)[i].index = k+j*numx;// ith particle is in cell # k+j*numx
+ (*bchain)[i].new_claimed_cells[0] = k+j*numx; // ith particle has just claimed cell number k+j*numx
+ (*bchain)[i].new_claimed_cells[1] = -1;// denotes end of claimed_cells list
+ // when we have finished claiming cells we call this function.
+ update_bchain2d( bchain, cells, i);
+ }//nump
+
+ claimed = 1;
+
+ while(claimed != 0){
+ claimed = 0 ;
+ for(i=0;i<nump;i++){
+ claim_cells2d( bchain, cells, pList, i);
+ claimed += (*bchain)[i].numclaimed;
+ update_bchain2d( bchain, cells, i);
+ }
+ }//while
+}
+void _DVCWeights_Calculate3D( void* dvcWeights, void* _swarm, Cell_LocalIndex lCell_I ) {
+ DVCWeights* self = (DVCWeights*) dvcWeights;
+ Swarm* swarm = (Swarm*) _swarm;
+ Particle_InCellIndex cParticleCount;
+ IntegrationPoint** particle;
+ static int visited = 0 ;
+ double dx,dy,dz,da;
+ static struct cell *cells;// the connected grid
+ struct particle *pList;// particle List
+ struct chain *bchain;//boundary chain
+ int nump,numx,numy,numz;
+ double BBXMIN = -1.0; // the ranges of the local coordinates of a FEM cell.
+ double BBXMAX = 1.0;
+ double BBYMIN = -1.0;
+ double BBYMAX = 1.0;
+ double BBZMIN = -1.0;
+ double BBZMAX = 1.0;
+ int i,k;
+
+ numx = self->resX;
+ numy = self->resY;
+ numz = self->resZ;
+
+ nump = cParticleCount = swarm->cellParticleCountTbl[lCell_I];
+
+ dx = (BBXMAX - BBXMIN)/numx;
+ dy = (BBYMAX - BBYMIN)/numy;
+ dz = (BBZMAX - BBZMIN)/numz;
+ da = dx*dy*dz;
+
+ // Construct the grid for the Voronoi cells only once.
+ // If we wanted to call this function again during a job with a different resolution
+ // then we should destroy the grid once we have looped through the whole mesh.
+ // I am assuming we are not going to do that for now.
+ // Easy to implement this anyway, if needed.
+ if(!visited){
+ /* The DVCWeights class should really be a class the next level up here */
+ /* We should be able to swap out the WeightsCalculator_CalculateAll instead of just setting
+ a pointer inside that function */
+ visited++;
+ construct_grid(&cells,numz,numy,numx,BBXMIN,BBYMIN,BBZMIN,BBXMAX,BBYMAX,BBZMAX);
+ }
+
+
+ // init the data structures
+ init_structs( &bchain, &pList, nump);
+ reset_grid(&cells,numz*numy*numx);
+
+ particle = (IntegrationPoint**)malloc(nump*sizeof(IntegrationPoint*));
+
+ // initialize the particle positions to be the local coordinates of the material swarm particles
+ // I am assuming the xi's (local coords) are precalculated somewhere and get reset based on material
+ // positions each time step.
+ for(i=0;i<nump;i++){
+
+ particle[i] = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, lCell_I, i );
+ pList[i].x = particle[i]->xi[0];
+ pList[i].y = particle[i]->xi[1];
+ pList[i].z = particle[i]->xi[2];
+
+ }
+ create_voronoi( &bchain, &pList, &cells, dx, dy, dz, nump, numx, numy, numz, BBXMIN, BBXMAX, BBYMIN, BBYMAX, BBZMIN, BBZMAX);
+ get_centroids( cells, pList,numz,numy,numx,nump,da);
+ // We are setting the integration points to be the centroids of the Voronoi regions here and
+ // the weight is the volume of each Voronoi region.
+ for(i=0;i<nump;i++){
+
+ particle[i]->xi[0] = pList[i].cx;
+ particle[i]->xi[1] = pList[i].cy;
+ particle[i]->xi[2] = pList[i].cz;
+ particle[i]->weight = pList[i].w;
+
+ }
+ for(k=0;k<nump;k++){
+ free(bchain[k].new_claimed_cells);
+ free(bchain[k].new_bound_cells);
+ }
+ free(particle);
+ free(bchain);
+ free(pList);
+
+}
+void _DVCWeights_Calculate2D( void* dvcWeights, void* _swarm, Cell_LocalIndex lCell_I ) {
+ DVCWeights* self = (DVCWeights*) dvcWeights;
+ Swarm* swarm = (Swarm*) _swarm;
+ Particle_InCellIndex cParticleCount;
+ IntegrationPoint** particle;
+ static int visited = 0 ;
+ double dx,dy,da;
+ static struct cell2d *cells;// the connected grid
+ struct particle2d *pList;// particle List
+ struct chain *bchain;//boundary chain
+ int nump,numx,numy;
+ double BBXMIN = -1.0; // the ranges of the local coordinates of a FEM cell.
+ double BBXMAX = 1.0;
+ double BBYMIN = -1.0;
+ double BBYMAX = 1.0;
+ int i,k;
+
+ numx = self->resX;
+ numy = self->resY;
+
+ nump = cParticleCount = swarm->cellParticleCountTbl[lCell_I];
+
+ dx = (BBXMAX - BBXMIN)/numx;
+ dy = (BBYMAX - BBYMIN)/numy;
+ da = dx*dy;
+
+ // Construct the grid for the Voronoi cells only once.
+ // If we wanted to call this function again during a job with a different resolution
+ // then we should destroy the grid once we have looped through the whole mesh.
+ // I am assuming we are not going to do that for now.
+ // Easy to implement this anyway, if needed.
+ if(!visited){
+ /* The DVCWeights class should really be a class the next level up here */
+ /* We should be able to swap out the WeightsCalculator_CalculateAll instead of just setting
+ a pointer inside that function */
+ visited++;
+ construct_grid2d(&cells,numy,numx,BBXMIN,BBYMIN,BBXMAX,BBYMAX);
+ }
+
+
+ // init the data structures
+ init_structs2d( &bchain, &pList, nump);
+ reset_grid2d(&cells,numy*numx);
+
+ particle = (IntegrationPoint**)malloc(nump*sizeof(IntegrationPoint*));
+
+ // initialize the particle positions to be the local coordinates of the material swarm particles
+ // I am assuming the xi's (local coords) are precalculated somewhere and get reset based on material
+ // positions each time step.
+ for(i=0;i<nump;i++){
+
+ particle[i] = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, lCell_I, i );
+ pList[i].x = particle[i]->xi[0];
+ pList[i].y = particle[i]->xi[1];
+
+ }
+ create_voronoi2d( &bchain, &pList, &cells, dx, dy, nump, numx, numy, BBXMIN, BBXMAX, BBYMIN, BBYMAX);
+ get_centroids2d( cells, pList,numy,numx,nump,da);
+ // We are setting the integration points to be the centroids of the Voronoi regions here and
+ // the weight is the volume of each Voronoi region.
+ for(i=0;i<nump;i++){
+
+ particle[i]->xi[0] = pList[i].cx;
+ particle[i]->xi[1] = pList[i].cy;
+ particle[i]->weight = pList[i].w;
+
+ }
+ for(k=0;k<nump;k++){
+ free(bchain[k].new_claimed_cells);
+ free(bchain[k].new_bound_cells);
+ }
+ free(particle);
+ free(bchain);
+ free(pList);
+
+}
+
+void _DVCWeights_Calculate( void* dvcWeights, void* _swarm, Cell_LocalIndex lCell_I ){
+ Swarm* swarm = (Swarm*) _swarm;
+ Dimension_Index dim = swarm->dim;
+
+ if(dim == 3){
+ _DVCWeights_Calculate3D( dvcWeights, _swarm, lCell_I);
+ }
+ else {
+ _DVCWeights_Calculate2D( dvcWeights, _swarm, lCell_I);
+ }
+}
+/*-------------------------------------------------------------------------------------------------------------------------
+** Public Functions
+*/
+
+
Added: long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.h 2006-10-31 21:33:19 UTC (rev 5159)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.h 2006-10-31 21:33:22 UTC (rev 5160)
@@ -0,0 +1,210 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+** Copyright (c) 2006, Monash Cluster Computing
+** All rights reserved.
+** Redistribution and use in source and binary forms, with or without modification,
+** are permitted provided that the following conditions are met:
+**
+** * Redistributions of source code must retain the above copyright notice,
+** this list of conditions and the following disclaimer.
+** * Redistributions in binary form must reproduce the above copyright
+** notice, this list of conditions and the following disclaimer in the
+** documentation and/or other materials provided with the distribution.
+** * Neither the name of the Monash University nor the names of its contributors
+** may be used to endorse or promote products derived from this software
+** without specific prior written permission.
+**
+** THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+** "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
+** THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+** PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS
+** BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+** CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+** SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+** OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+**
+**
+** Contact:
+*% Louis Moresi - Louis.Moresi at sci.monash.edu.au
+*%
+** Author:
+** Mirko Velic - Mirko.Velic at sci.monash.edu.au
+**
+** Assumptions:
+** I am assuming that the xi's (local coords) on the IntegrationPoint particles
+** are precalculated somewhere and get reset based on material PIC positions each time step.
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+
+/****************************************************************************************************************
+
+ The algorithm here-in computes a discrete voronoi diagram per FEM cell given a set of local
+ particle positions, in 3D and 2D. The volumes of the Voronoi regions are used as integration weights for
+ the integration point swarm and the integration points are the centroids of the same volumes.
+
+ For a description of this algorithm, see the article by Velic et.al.
+ "A Fast Robust Algorithm for computing Discrete Voronoi Diagrams in N-dimensions"
+
+
+
+*****************************************************************************************************************/
+
+#ifndef __Underworld_Utils_DVCWeightsClass_h__
+#define __Underworld_Utils_DVCWeightsClass_h__
+
+ /* Textual name of this class */
+ extern const Type DVCWeights_Type;
+
+ /* DVCWeights information */
+ #define __DVCWeights \
+ /* General info */ \
+ __WeightsCalculator \
+ /* Virtual Info */\
+ /* Parameters that are passed in */ \
+ int resX; \
+ int resY; \
+ int resZ;
+
+ struct DVCWeights { __DVCWeights };
+
+#define INC 50
+
+struct cell{
+ int p;//particle index number
+ int index;
+ int N;
+ int S;
+ int E;
+ int W;
+ int U;
+ int D;
+ double x;
+ double y;
+ double z;
+ int done;
+};
+struct cell2d{
+ int p;//particle index number
+ int index;
+ int N;
+ int S;
+ int E;
+ int W;
+ double x;
+ double y;
+ int done;
+};
+struct chain{
+ int p;//particle index number
+ int index;//cell number in grid
+ int sizeofboundary; // number of cells on boundary so far
+ int numclaimed;
+ int totalclaimed;
+ int new_bound_cells_malloced;
+ int new_claimed_cells_malloced;
+ int *new_bound_cells;
+ int *new_claimed_cells;
+ int done;
+};
+struct particle{
+ double x;
+ double y;
+ double z;
+ double cx;
+ double cy;
+ double cz;
+ double w;
+ int index;
+};
+struct particle2d{
+ double x;
+ double y;
+ double cx;
+ double cy;
+ double w;
+ int index;
+};
+// function prototypes
+void get_centroids( struct cell *cells,struct particle *pList,
+ int n,int m, int l,int nump,double vol);
+void claim_cells(struct chain **bchain,struct cell **cells,struct particle **pList,int p_i);
+void update_bchain(struct chain **bchain,struct cell **cells,int p_i);
+void reset_grid(struct cell **cells, int n );
+void reset_grid2d(struct cell2d **cells, int n );
+double distance(double x0, double y0, double z0, double x1, double y1, double z1);
+void construct_grid(struct cell **cell_list, int n,int m, int l,
+ double x0,double y0,double z0,double x1,double y1,double z1);
+void init_structs( struct chain **bchain, struct particle **pList, int nump);
+void create_voronoi( struct chain **bchain, struct particle **pList, struct cell **cells,
+ double dx, double dy, double dz,
+ int nump,
+ int numx, int numy, int numz,
+ double BBXMIN, double BBXMAX,
+ double BBYMIN, double BBYMAX,
+ double BBZMIN, double BBZMAX);
+
+void get_centroids2d( struct cell2d *cells, struct particle2d *pList,
+ int n, int m, int nump,double vol);
+void claim_cells2d(struct chain **bchain,struct cell2d **cells,struct particle2d **pList,int p_i);
+void update_bchain2d(struct chain **bchain,struct cell2d **cells,int p_i);
+
+double distance2d(double x0, double y0, double x1, double y1);
+void construct_grid2d(struct cell2d **cell_list, int m, int l,
+ double x0,double y0,double x1,double y1);
+void init_structs2d( struct chain **bchain, struct particle2d **pList, int nump);
+void create_voronoi2d( struct chain **bchain, struct particle2d **pList, struct cell2d **cells,
+ double dx, double dy,
+ int nump,
+ int numx, int numy,
+ double BBXMIN, double BBXMAX,
+ double BBYMIN, double BBYMAX);
+
+void _DVCWeights_Calculate2D( void* dvcWeights, void* _swarm, Cell_LocalIndex lCell_I );
+void _DVCWeights_Calculate3D( void* dvcWeights, void* _swarm, Cell_LocalIndex lCell_I );
+// end prototypes
+
+ /*---------------------------------------------------------------------------------------------------------------------
+ ** Constructors
+ */
+ DVCWeights* DVCWeights_New( Name name, Dimension_Index dim ) ;
+ DVCWeights* _DVCWeights_New(
+ SizeT _sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ WeightsCalculator_CalculateFunction* _calculate,
+ Name name );
+
+ void _DVCWeights_Init( void* dvcWeights, int *res ) ;
+ void DVCWeights_InitAll( void* dvcWeights, Dimension_Index dim ) ;
+
+
+ /* Stg_Class_Delete DVCWeights implementation */
+ void _DVCWeights_Delete( void* dvcWeights );
+ void _DVCWeights_Print( void* dvcWeights, Stream* stream );
+ #define DVCWeights_Copy( self ) \
+ (DVCWeights*) Stg_Class_Copy( self, NULL, False, NULL, NULL )
+ #define DVCWeights_DeepCopy( self ) \
+ (DVCWeights*) Stg_Class_Copy( self, NULL, True, NULL, NULL )
+ void* _DVCWeights_Copy( void* dvcWeights, void* dest, Bool deep, Name nameExt, PtrMap* ptrMap );
+
+ void* _DVCWeights_DefaultNew( Name name ) ;
+ void _DVCWeights_Construct( void* dvcWeights, Stg_ComponentFactory* cf ) ;
+ void _DVCWeights_Build( void* dvcWeights, void* data ) ;
+ void _DVCWeights_Initialise( void* dvcWeights, void* data ) ;
+ void _DVCWeights_Execute( void* dvcWeights, void* data );
+ void _DVCWeights_Destroy( void* dvcWeights, void* data ) ;
+
+ void _DVCWeights_Calculate( void* dvcWeights, void* _swarm, Cell_LocalIndex lCell_I ) ;
+
+
+#endif
Added: long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.meta
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.meta 2006-10-31 21:33:19 UTC (rev 5159)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/DVCWeights.meta 2006-10-31 21:33:22 UTC (rev 5160)
@@ -0,0 +1,54 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+<param name="Name">DVCWeights</param>
+<param name="Organisation">MCC</param>
+<param name="Project">Underworld</param>
+<param name="Location">./Underworld/Utils/src/</param>
+<param name="Project Web">http://mcc.monash.edu.au/Underworld</param>
+<param name="Copyright">Copyright (C) 2006 Monash Cluster Computing.</param>
+<param name="License">http://www.opensource.org/licenses/bsd-license.php</param>
+<param name="Parent">WeightsCalculator</param>
+<param name="Description">Calculates PIC integration weights in 3D and 2D using a Voronoi diagram. The integration points are set to the centroid of a Voronoi cell and
+ the weight is the volume of the cell.</param>
+
+
+<list name="Params">
+
+</list>
+
+<list name="Dependencies">
+
+</list>
+<!-- Example Usage -->
+<struct name="weights">
+ <param name="Type">DVCWeights</param>
+ <param name="resolution">10</param>
+</struct>
+<!-- or one can call this with -->
+<struct name="weights">
+ <param name="Type">DVCWeights</param>
+ <param name="resolutionX">10</param>
+ <param name="resolutionY">10</param>
+ <param name="resolutionZ">10</param>
+</struct>
+
+<struct name="materialPoints">
+ <param name="Type">IntegrationPointsSwarm</param>
+ <param name="CellLayout">elementCellLayout</param>
+ <param name="ParticleLayout">localLayout</param>
+ <param name="FiniteElement_Mesh">mesh-linear</param>
+ <param name="WeightsCalculator">weights</param>
+ <param name="TimeIntegrator">timeIntegrator</param>
+ <param name="IntegrationPointMapper">mapper</param>
+</struct>
+<struct name="mapper">
+ <param name="Type">CoincidentMapper</param>
+ <param name="IntegrationPointsSwarm">materialPoints</param>
+ <param name="MaterialPointsSwarm">materialSwarm</param>
+</struct>
+
+<!-- See MaterialPointsSwarm.xml in Underworld/InputFiles/PIC_Components for an example -->
+<!-- where the X Y Z resolutions may be all different -->
+</StGermainData>
Modified: long/3D/Gale/trunk/src/Underworld/Utils/src/Init.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/Init.c 2006-10-31 21:33:19 UTC (rev 5159)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/Init.c 2006-10-31 21:33:22 UTC (rev 5160)
@@ -46,7 +46,11 @@
#include <mpi.h>
#include <StGermain/StGermain.h>
#include <StgFEM/StgFEM.h>
+
#include <PICellerator/PICellerator.h>
+#include <PICellerator/Weights/Weights.h>
+#include <PICellerator/Weights/WeightsCalculator.h>
+
#include <Underworld/Rheology/Rheology.h>
#include "Utils.h"
@@ -63,12 +67,14 @@
Stg_ComponentRegister_Add( componentRegister, RadiogenicHeatingTerm_Type, "0", _RadiogenicHeatingTerm_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, StressField_Type , "0", _StressField_DefaultNew );
Stg_ComponentRegister_Add( componentRegister, ViscosityField_Type , "0", _ViscosityField_DefaultNew );
-
+ Stg_ComponentRegister_Add( componentRegister, DVCWeights_Type, "0", _DVCWeights_DefaultNew );
+
RegisterParent( UnderworldContext_Type, PICelleratorContext_Type );
RegisterParent( PressureTemperatureOutput_Type, SwarmOutput_Type );
RegisterParent( RadiogenicHeatingTerm_Type, ForceTerm_Type );
RegisterParent( StressField_Type, ParticleFeVariable_Type );
RegisterParent( ViscosityField_Type, ParticleFeVariable_Type );
+ RegisterParent( DVCWeights_Type, WeightsCalculator_Type );
return True;
}
Modified: long/3D/Gale/trunk/src/Underworld/Utils/src/Utils.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/Utils.h 2006-10-31 21:33:19 UTC (rev 5159)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/Utils.h 2006-10-31 21:33:22 UTC (rev 5160)
@@ -50,6 +50,7 @@
#include "Context.h"
#include "StressField.h"
#include "ViscosityField.h"
+ #include "DVCWeights.h"
#include "PressureTemperatureOutput.h"
#include "RadiogenicHeatingTerm.h"
#include "Init.h"
Modified: long/3D/Gale/trunk/src/Underworld/Utils/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/types.h 2006-10-31 21:33:19 UTC (rev 5159)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/types.h 2006-10-31 21:33:22 UTC (rev 5160)
@@ -51,5 +51,6 @@
typedef struct RadiogenicHeatingTerm RadiogenicHeatingTerm;
typedef struct StressField StressField;
typedef struct ViscosityField ViscosityField;
+ typedef struct DVCWeights DVCWeights;
#endif
More information about the cig-commits
mailing list