[cig-commits] commit: Add BuoyancyDampingTerm
Mercurial
hg at geodynamics.org
Mon Apr 30 13:50:57 PDT 2012
changeset: 449:c1e2c18c39f8
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Mon Apr 30 13:50:43 2012 -0700
files: Utils/src/BuoyancyDampingTerm.cxx Utils/src/BuoyancyDampingTerm.h Utils/src/BuoyancyDampingTerm.meta Utils/src/BuoyancyForceTerm.cxx Utils/src/BuoyancyForceTerm.h Utils/src/Init.cxx Utils/src/Utils.h Utils/src/types.h
description:
Add BuoyancyDampingTerm
diff -r 30411e58fae6 -r c1e2c18c39f8 Utils/src/BuoyancyDampingTerm.cxx
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Utils/src/BuoyancyDampingTerm.cxx Mon Apr 30 13:50:43 2012 -0700
@@ -0,0 +1,295 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
+** Melbourne, 3053, Australia.
+**
+** Primary Contributing Organisations:
+** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
+** Australian Computational Earth Systems Simulator - http://www.access.edu.au
+** Monash Cluster Computing - http://www.mcc.monash.edu.au
+** Computational Infrastructure for Geodynamics - http://www.geodynamics.org
+**
+** Contributors:
+** Patrick D. Sunter, Software Engineer, VPAC. (pds at vpac.org)
+** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
+** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
+** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+** $Id: BuoyancyDampingTerm.c 964 2007-10-11 08:03:06Z SteveQuenette $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <assert.h>
+#include <string.h>
+#include <mpi.h>
+#include <StGermain/StGermain.h>
+#include <StgDomain/StgDomain.h>
+#include <StgFEM/StgFEM.h>
+#include <PICellerator/PopulationControl/PopulationControl.h>
+#include <PICellerator/Weights/Weights.h>
+#include <PICellerator/MaterialPoints/MaterialPoints.h>
+#include "types.h"
+#include "BuoyancyForceTerm.h"
+#include "BuoyancyDampingTerm.h"
+#include "MaterialSwarmVariable.h"
+
+
+/* Textual name of this class */
+const Type BuoyancyDampingTerm_Type = "BuoyancyDampingTerm";
+
+
+BuoyancyDampingTerm*
+BuoyancyDampingTerm_New(Name name,
+ FiniteElementContext* context,
+ StiffnessMatrix* stiffnessMatrix,
+ Swarm* integrationSwarm,
+ BuoyancyForceTerm* buoyancyForceTerm)
+{
+ BuoyancyDampingTerm* self =
+ (BuoyancyDampingTerm*) _BuoyancyDampingTerm_DefaultNew( name );
+
+ self->isConstructed = True;
+ _StiffnessMatrixTerm_Init(self,context,stiffnessMatrix,integrationSwarm,NULL);
+
+ _BuoyancyDampingTerm_Init(self, buoyancyForceTerm);
+
+ return self;
+}
+
+/* Creation implementation / Virtual constructor */
+BuoyancyDampingTerm* _BuoyancyDampingTerm_New( BUOYANCYDAMPINGTERM_DEFARGS )
+{
+ BuoyancyDampingTerm* self;
+
+ /* Allocate memory */
+ assert( _sizeOfSelf >= sizeof(BuoyancyDampingTerm) );
+ self=(BuoyancyDampingTerm*)
+ _StiffnessMatrixTerm_New( STIFFNESSMATRIXTERM_PASSARGS );
+
+ /* Virtual info */
+
+ return self;
+}
+
+void _BuoyancyDampingTerm_Init(void* matrixTerm,
+ BuoyancyForceTerm* buoyancyForceTerm)
+{
+ BuoyancyDampingTerm* self = (BuoyancyDampingTerm*)matrixTerm;
+ self->buoyancyForceTerm=buoyancyForceTerm;
+}
+
+void _BuoyancyDampingTerm_Delete( void* matrixTerm ) {
+ BuoyancyDampingTerm* self = (BuoyancyDampingTerm*)matrixTerm;
+ _StiffnessMatrixTerm_Delete( self );
+}
+
+void _BuoyancyDampingTerm_Print( void* matrixTerm, Stream* stream ) {
+ BuoyancyDampingTerm* self = (BuoyancyDampingTerm*)matrixTerm;
+
+ _StiffnessMatrixTerm_Print( self, stream );
+
+ /* General info */
+}
+
+void* _BuoyancyDampingTerm_DefaultNew(Name name)
+{
+ /* Variables set in this function */
+ SizeT _sizeOfSelf = sizeof(BuoyancyDampingTerm);
+ Type type = BuoyancyDampingTerm_Type;
+ Stg_Class_DeleteFunction* _delete= _BuoyancyDampingTerm_Delete;
+ Stg_Class_PrintFunction* _print= _BuoyancyDampingTerm_Print;
+ Stg_Class_CopyFunction* _copy= NULL;
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor=
+ _BuoyancyDampingTerm_DefaultNew;
+ Stg_Component_ConstructFunction* _construct=_BuoyancyDampingTerm_AssignFromXML;
+ Stg_Component_BuildFunction* _build=_BuoyancyDampingTerm_Build;
+ Stg_Component_InitialiseFunction* _initialise=_BuoyancyDampingTerm_Initialise;
+ Stg_Component_ExecuteFunction* _execute=_BuoyancyDampingTerm_Execute;
+ Stg_Component_DestroyFunction* _destroy=_BuoyancyDampingTerm_Destroy;
+ StiffnessMatrixTerm_AssembleElementFunction* _assembleElement=
+ _BuoyancyDampingTerm_AssembleElement;
+
+ /* Variables that are set to ZERO are variables that will be
+ set either by the current _New function or another parent
+ _New function further up the hierachy */
+ AllocationType nameAllocationType = NON_GLOBAL /* default
+ value
+ NON_GLOBAL */;
+
+ return (void*)_BuoyancyDampingTerm_New( BUOYANCYDAMPINGTERM_PASSARGS );
+}
+
+void _BuoyancyDampingTerm_AssignFromXML(void* matrixTerm,
+ Stg_ComponentFactory* cf, void* data)
+{
+ BuoyancyDampingTerm* self=(BuoyancyDampingTerm*)matrixTerm;
+ /* Construct Parent */
+ _StiffnessMatrixTerm_AssignFromXML( self, cf, data );
+
+ BuoyancyForceTerm *force=(BuoyancyForceTerm*)
+ Stg_ComponentFactory_ConstructByKey(cf,self->name,(Dictionary_Entry_Key)"BuoyancyForceTerm",
+ BuoyancyForceTerm,False,data);
+ _BuoyancyDampingTerm_Init(self, force);
+}
+
+void _BuoyancyDampingTerm_Build( void* matrixTerm, void* data )
+{
+ BuoyancyDampingTerm* self = (BuoyancyDampingTerm*)matrixTerm;
+ _StiffnessMatrixTerm_Build( self, data );
+ Stg_Component_Build( self->buoyancyForceTerm, data, False);
+}
+
+void _BuoyancyDampingTerm_Initialise( void* matrixTerm, void* data )
+{
+ BuoyancyDampingTerm* self = (BuoyancyDampingTerm*)matrixTerm;
+
+ _StiffnessMatrixTerm_Initialise( self, data );
+ Stg_Component_Initialise( self->buoyancyForceTerm, data, False );
+}
+
+void _BuoyancyDampingTerm_Execute( void* matrixTerm, void* data )
+{
+ _StiffnessMatrixTerm_Execute( matrixTerm, data );
+}
+
+void _BuoyancyDampingTerm_Destroy( void* matrixTerm, void* data )
+{
+ _StiffnessMatrixTerm_Destroy( matrixTerm, data );
+}
+
+
+void _BuoyancyDampingTerm_AssembleElement(void* matrixTerm,
+ StiffnessMatrix* stiffnessMatrix,
+ Element_LocalIndex lElement_I,
+ SystemLinearEquations* sle,
+ FiniteElementContext* context,
+ double** elStiffMat)
+{
+ BuoyancyDampingTerm* self=Stg_CheckType(matrixTerm,BuoyancyDampingTerm);
+ if(!self->buoyancyForceTerm || !self->buoyancyForceTerm->damping)
+ return;
+
+ FeVariable* variable1=stiffnessMatrix->rowVariable;
+ Dimension_Index dim=stiffnessMatrix->dim;
+ Dof_Index nodeDofCount=dim;
+
+ /* Set the element type */
+ ElementType* elementType=FeMesh_GetElementType(variable1->feMesh,lElement_I);
+
+ double density, alpha;
+
+ BuoyancyForceTerm_average_density_alpha
+ (self->buoyancyForceTerm,
+ (IntegrationPointsSwarm*)self->buoyancyForceTerm->integrationSwarm,
+ variable1->feMesh,lElement_I,&density,&alpha);
+
+ double **jac=Memory_Alloc_2DArray(double,dim,dim,(Name)"Temporary Jacobian");
+
+ /* Gauss-Legendre integration */
+ const double legendre_points[]={-sqrt(3/5.0),0,sqrt(3/5.0)};
+ const double weights[]={5/9.0, 8/9.0, 5/9.0};
+
+ for(unsigned int local_norm=0;local_norm<dim; ++local_norm)
+ {
+ const int tangent((local_norm+1)%dim), tangent2((tangent+1)%dim);
+
+ for(int sgn=-1;sgn<=1;sgn+=2)
+ {
+ const unsigned int face((local_norm==0 ? 1 :
+ (local_norm==1 ? 0 : 2))*2+(sgn+1)/2);
+ for(int i=0; i<3; ++i)
+ for(int j=0;j<(dim==2 ? 1 : 3); ++j)
+ {
+ double xi[dim];
+ xi[local_norm]=sgn;
+ xi[tangent]=legendre_points[i];
+ if(dim==3)
+ {
+ xi[tangent2]=legendre_points[j];
+ }
+ const int num_nodes(9);
+ double Ni[num_nodes];
+ ElementType_EvaluateShapeFunctionsAt(elementType,xi,Ni);
+ ElementType_Jacobian(elementType,variable1->feMesh,
+ lElement_I,xi,dim,jac,NULL);
+
+ double gravity[dim];
+ BuoyancyForceTerm_CalcGravity(self->buoyancyForceTerm,
+ variable1->feMesh,lElement_I,
+ xi,gravity);
+ double temperature(0);
+ if(alpha!=0 && self->buoyancyForceTerm->temperatureField)
+ FeVariable_InterpolateFromMeshLocalCoord
+ (self->buoyancyForceTerm->temperatureField,
+ variable1->feMesh, lElement_I, xi, &temperature);
+
+ double rho(density*(1.0-alpha*temperature));
+ for(unsigned int d=0;d<dim;++d)
+ {
+ double damping(-gravity[d]*rho*context->dt);
+ double geometric_factor;
+
+ /* For a top boundary, consider y=a*x. So
+ dx/dxi>0, and dx/deta=1/a. If a>0, Fy<0, so
+ need to change the sign with the d==local_norm
+ term. */
+
+ if(dim==2)
+ {
+ geometric_factor=jac[tangent][(d+1)%dim]
+ *sgn*(d==local_norm ? 1 : -1);
+ }
+ else
+ {
+ geometric_factor=
+ (jac[tangent][(d+1)%dim] * jac[tangent2][(d+2)%dim]
+ -jac[tangent][(d+2)%dim]*jac[tangent2][(d+1)%dim])
+ *sgn*(d==local_norm ? 1 : -1);
+ }
+ /* Now loop over the nodes of the element */
+ const int num_face_nodes(dim==2 ? (num_nodes==4 ? 2 : 3)
+ : (num_nodes==8 ? 4 : 9));
+ for(int el_0=0;el_0<num_face_nodes;++el_0)
+ {
+ const int node_0(elementType->faceNodes[face][el_0]);
+ for(int el_1=0;el_1<num_face_nodes;++el_1)
+ {
+ const int node_1
+ (elementType->faceNodes[face][el_1]);
+
+ double total=damping*weights[i]*geometric_factor
+ *Ni[node_0]*Ni[node_1];
+ if(dim==3)
+ {
+ total*=weights[j];
+ }
+ elStiffMat[node_0*nodeDofCount+d]
+ [node_1*nodeDofCount+d] +=total;
+ }
+ }
+ }
+ }
+ }
+ }
+ Memory_Free( jac );
+}
diff -r 30411e58fae6 -r c1e2c18c39f8 Utils/src/BuoyancyDampingTerm.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Utils/src/BuoyancyDampingTerm.h Mon Apr 30 13:50:43 2012 -0700
@@ -0,0 +1,109 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
+** Melbourne, 3053, Australia.
+**
+** Primary Contributing Organisations:
+** Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
+** Australian Computational Earth Systems Simulator - http://www.access.edu.au
+** Monash Cluster Computing - http://www.mcc.monash.edu.au
+** Computational Infrastructure for Geodynamics - http://www.geodynamics.org
+**
+** Contributors:
+** Patrick D. Sunter, Software Engineer, VPAC. (pds at vpac.org)
+** Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
+** Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
+** Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+*/
+
+
+#ifndef __PICellerator_Utils_BuoyancyDampingTerm_h__
+#define __PICellerator_Utils_BuoyancyDampingTerm_h__
+
+/** Textual name of this class */
+extern const Type BuoyancyDampingTerm_Type;
+
+ /** BuoyancyDampingTerm class contents */
+ #define __BuoyancyDampingTerm \
+ /* General info */ \
+ __StiffnessMatrixTerm \
+ \
+ /* Virtual info */ \
+ \
+ /* BuoyancyDampingTerm info */ \
+ \
+ BuoyancyForceTerm* buoyancyForceTerm;
+
+ struct BuoyancyDampingTerm { __BuoyancyDampingTerm };
+
+ BuoyancyDampingTerm* BuoyancyDampingTerm_New(
+ Name name,
+ FiniteElementContext* context,
+ StiffnessMatrix* stiffnessMatrix,
+ Swarm* integrationSwarm,
+ BuoyancyForceTerm* buoyancyForceTerm);
+
+
+ #ifndef ZERO
+ #define ZERO 0
+ #endif
+
+ #define BUOYANCYDAMPINGTERM_DEFARGS \
+ STIFFNESSMATRIXTERM_DEFARGS
+
+ #define BUOYANCYDAMPINGTERM_PASSARGS \
+ STIFFNESSMATRIXTERM_PASSARGS
+
+ BuoyancyDampingTerm* _BuoyancyDampingTerm_New( BUOYANCYDAMPINGTERM_DEFARGS );
+
+ void _BuoyancyDampingTerm_Init(void* matrixTerm,
+ BuoyancyForceTerm* buoyancyForceTerm);
+
+ void _BuoyancyDampingTerm_Delete( void* matrixTerm );
+
+ void _BuoyancyDampingTerm_Print( void* matrixTerm, Stream* stream );
+
+ void* _BuoyancyDampingTerm_DefaultNew( Name name );
+
+ void _BuoyancyDampingTerm_AssignFromXML( void* matrixTerm,
+ Stg_ComponentFactory* cf,
+ void* data );
+
+ void _BuoyancyDampingTerm_Build( void* matrixTerm, void* data );
+
+ void _BuoyancyDampingTerm_Initialise( void* matrixTerm, void* data );
+
+ void _BuoyancyDampingTerm_Execute( void* matrixTerm, void* data );
+
+ void _BuoyancyDampingTerm_Destroy( void* matrixTerm, void* data );
+
+ void _BuoyancyDampingTerm_AssembleElement(void* matrixTerm,
+ StiffnessMatrix* stiffnessMatrix,
+ Element_LocalIndex lElement_I,
+ SystemLinearEquations* sle,
+ FiniteElementContext* context,
+ double** elStiffMat ) ;
+
+#endif
+
diff -r 30411e58fae6 -r c1e2c18c39f8 Utils/src/BuoyancyDampingTerm.meta
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/Utils/src/BuoyancyDampingTerm.meta Mon Apr 30 13:50:43 2012 -0700
@@ -0,0 +1,31 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+<param name="Name">BuoyancyDampingTerm</param>
+<param name="Author">...</param>
+<param name="Organisation">VPAC</param>
+<param name="Project">StgFEM</param>
+<param name="Location">./PICellerator/Utils/src/</param>
+<param name="Project Web">http://www.stgermainproject.org/StgFEM.html</param>
+<param name="Copyright">Copyright (C) 2004-2005 VPAC.</param>
+<param name="License">The Gnu Lesser General Public License v2.1 - http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html</param>
+<param name="Parent">StiffnessMatrixTerm</param>
+<param name="Reference">...</param>
+<param name="Summary">...</param>
+<param name="Description">...</param>
+
+<!--Now the interesting stuff-->
+
+
+<list name="Params">
+
+</list>
+
+<list name="Dependencies">
+
+</list>
+<!-- Add an exmaple XML if possible -->
+<param name="Example">...</param>
+
+</StGermainData>
diff -r 30411e58fae6 -r c1e2c18c39f8 Utils/src/BuoyancyForceTerm.cxx
--- a/Utils/src/BuoyancyForceTerm.cxx Sat Apr 07 15:10:18 2012 -0700
+++ b/Utils/src/BuoyancyForceTerm.cxx Mon Apr 30 13:50:43 2012 -0700
@@ -64,7 +64,7 @@
#include <stddef.h>
#include <sstream>
-static std::vector<std::string> densities, alphas;
+static std::vector<std::string> equations[2];
/* Textual name of this class */
const Type BuoyancyForceTerm_Type = "BuoyancyForceTerm";
@@ -77,7 +77,9 @@ BuoyancyForceTerm* BuoyancyForceTerm_New
FeVariable* temperatureField,
FeVariable* pressureField,
double gravity,
+ double* gHat,
Bool adjust,
+ Bool damping,
Materials_Register* materials_Register,
HydrostaticTerm* hydrostaticTerm )
{
@@ -85,8 +87,9 @@ BuoyancyForceTerm* BuoyancyForceTerm_New
self->isConstructed = True;
_ForceTerm_Init( self, context, forceVector, integrationSwarm, NULL );
- _BuoyancyForceTerm_Init( self, temperatureField, pressureField, gravity, adjust, materials_Register,
- hydrostaticTerm );
+ _BuoyancyForceTerm_Init(self,temperatureField,pressureField,gravity,
+ gHat,adjust,damping,materials_Register,
+ hydrostaticTerm);
return self;
}
@@ -116,7 +119,9 @@ void _BuoyancyForceTerm_Init(
FeVariable* temperatureField,
FeVariable* pressureField,
double gravity,
+ double* gHat,
Bool adjust,
+ Bool damping,
Materials_Register* materials_Register,
HydrostaticTerm* hydrostaticTerm )
{
@@ -125,8 +130,9 @@ void _BuoyancyForceTerm_Init(
self->temperatureField = temperatureField;
self->pressureField = pressureField;
self->gravity = gravity;
- self->gHat = NULL;
+ self->gHat = gHat;
self->adjust = adjust;
+ self->damping = damping;
self->materials_Register = materials_Register;
self->hydrostaticTerm = hydrostaticTerm;
}
@@ -170,193 +176,239 @@ void* _BuoyancyForceTerm_DefaultNew( Nam
return (void*)_BuoyancyForceTerm_New( BUOYANCYFORCETERM_PASSARGS );
}
-void _BuoyancyForceTerm_AssignFromXML( void* forceTerm, Stg_ComponentFactory* cf, void* data ) {
- BuoyancyForceTerm* self = (BuoyancyForceTerm*)forceTerm;
- Dictionary* dict;
- Dictionary_Entry_Value* tmp;
- char* rootKey;
- FeVariable* temperatureField;
- FeVariable* pressureField;
- double gravity;
- Bool adjust;
- Materials_Register* materials_Register;
- unsigned nDims;
- Dictionary_Entry_Value* direcList;
- double* direc;
- unsigned d_i;
- PICelleratorContext* context;
- HydrostaticTerm* hydrostaticTerm;
+void _BuoyancyForceTerm_AssignFromXML(void* forceTerm,Stg_ComponentFactory* cf,
+ void* data)
+{
+ BuoyancyForceTerm* self = (BuoyancyForceTerm*)forceTerm;
+ Dictionary* dict;
+ Dictionary_Entry_Value* tmp;
+ char* rootKey;
+ FeVariable* temperatureField;
+ FeVariable* pressureField;
+ double gravity;
+ Bool adjust;
+ Bool damping;
+ Materials_Register* materials_Register;
+ unsigned nDims;
+ Dictionary_Entry_Value* direcList;
+ double* gHat;
+ unsigned d_i;
+ PICelleratorContext* context;
+ HydrostaticTerm* hydrostaticTerm;
- /* Construct Parent */
- _ForceTerm_AssignFromXML( self, cf, data );
+ /* Construct Parent */
+ _ForceTerm_AssignFromXML( self, cf, data );
- dict = Dictionary_Entry_Value_AsDictionary( Dictionary_Get( cf->componentDict, (Dictionary_Entry_Key)self->name ) );
- temperatureField = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"TemperatureField", FeVariable, False, data ) ;
- pressureField = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"PressureField", FeVariable, False, data ) ;
- gravity = Stg_ComponentFactory_GetDouble( cf, self->name, (Dictionary_Entry_Key)"gravity", 0.0 );
- adjust = Stg_ComponentFactory_GetBool( cf, self->name, (Dictionary_Entry_Key)"adjust", False );
+ dict=Dictionary_Entry_Value_AsDictionary(Dictionary_Get
+ (cf->componentDict,
+ (Dictionary_Entry_Key)self->name));
+
+ temperatureField=
+ Stg_ComponentFactory_ConstructByKey(cf,self->name,
+ (Dictionary_Entry_Key)"TemperatureField",
+ FeVariable, False, data ) ;
+ pressureField=
+ Stg_ComponentFactory_ConstructByKey(cf,self->name,
+ (Dictionary_Entry_Key)"PressureField",
+ FeVariable,False,data);
+ gravity=
+ Stg_ComponentFactory_GetDouble(cf,self->name,
+ (Dictionary_Entry_Key)"gravity",0.0);
+ adjust=Stg_ComponentFactory_GetBool(cf,self->name,
+ (Dictionary_Entry_Key)"adjust",False);
+ damping=Stg_ComponentFactory_GetBool(cf,self->name,
+ (Dictionary_Entry_Key)"damping",True);
- direcList = Dictionary_Get( dict, (Dictionary_Entry_Key)"gravityDirection" );
+ direcList=Dictionary_Get(dict,(Dictionary_Entry_Key)"gravityDirection");
- if( direcList ) {
- nDims = Dictionary_Entry_Value_GetCount( direcList );
- direc = AllocArray( double, nDims );
+ if(direcList)
+ {
+ nDims = Dictionary_Entry_Value_GetCount( direcList );
+ gHat = AllocArray( double, nDims );
- for( d_i = 0; d_i < nDims; d_i++ ) {
- tmp = Dictionary_Entry_Value_GetElement( direcList, d_i );
- rootKey = Dictionary_Entry_Value_AsString( tmp );
+ for( d_i = 0; d_i < nDims; d_i++ )
+ {
+ tmp = Dictionary_Entry_Value_GetElement( direcList, d_i );
+ rootKey = Dictionary_Entry_Value_AsString( tmp );
- if( !Stg_StringIsNumeric( (char *)rootKey ) )
- tmp = Dictionary_Get( cf->rootDict, (Dictionary_Entry_Key)rootKey );
- direc[d_i] = Dictionary_Entry_Value_AsDouble( tmp );
- }
- if( nDims == 2 )
- Vec_Norm2D( direc, direc );
- else
- Vec_Norm3D( direc, direc );
- }
- else
- direc = NULL;
- self->gHat = direc;
+ if( !Stg_StringIsNumeric( (char *)rootKey ) )
+ tmp = Dictionary_Get( cf->rootDict, (Dictionary_Entry_Key)rootKey );
+ gHat[d_i] = Dictionary_Entry_Value_AsDouble( tmp );
+ }
+ if( nDims == 2 )
+ Vec_Norm2D( gHat, gHat );
+ else
+ Vec_Norm3D( gHat, gHat );
+ }
+ else
+ {
+ gHat = AllocArray( double, 3 );
+ gHat[0]=0;
+ gHat[1]=-1;
+ gHat[2]=0;
+ }
- context = (PICelleratorContext*)self->context;
- assert( Stg_CheckType( context, PICelleratorContext ) );
- materials_Register = context->materials_Register;
- assert( materials_Register );
+ context = (PICelleratorContext*)self->context;
+ assert( Stg_CheckType( context, PICelleratorContext ) );
+ materials_Register = context->materials_Register;
+ assert( materials_Register );
- hydrostaticTerm = Stg_ComponentFactory_ConstructByKey( cf, self->name, "HydrostaticTerm", HydrostaticTerm, False, data ) ;
+ hydrostaticTerm = Stg_ComponentFactory_ConstructByKey(cf, self->name,
+ "HydrostaticTerm",
+ HydrostaticTerm, False,
+ data);
- _BuoyancyForceTerm_Init( self, temperatureField, pressureField,
- gravity, adjust,
- materials_Register, hydrostaticTerm );
+ _BuoyancyForceTerm_Init(self,temperatureField,pressureField,gravity,gHat,
+ adjust,damping,materials_Register,hydrostaticTerm);
}
-void _BuoyancyForceTerm_Build( void* forceTerm, void* data ) {
- BuoyancyForceTerm* self = (BuoyancyForceTerm*)forceTerm;
- BuoyancyForceTerm_MaterialExt* materialExt;
- Material_Index material_I;
- Material* material;
- Materials_Register* materials_Register = self->materials_Register;
- IntegrationPointsSwarm* swarm = (IntegrationPointsSwarm*)self->integrationSwarm;
- MaterialPointsSwarm** materialSwarms;
- Index materialSwarm_I;
- char* name;
- Stg_ComponentFactory* cf;
+void _BuoyancyForceTerm_Build(void* forceTerm, void* data)
+{
+ BuoyancyForceTerm* self = (BuoyancyForceTerm*)forceTerm;
+ BuoyancyForceTerm_MaterialExt* materialExt;
+ Material_Index material_I;
+ Material* material;
+ Materials_Register* materials_Register = self->materials_Register;
+ IntegrationPointsSwarm* swarm=(IntegrationPointsSwarm*)self->integrationSwarm;
+ MaterialPointsSwarm** materialSwarms;
+ Index materialSwarm_I;
+ char* name;
+ Stg_ComponentFactory* cf;
- cf = self->context->CF;
+ cf = self->context->CF;
- _ForceTerm_Build( self, data );
+ _ForceTerm_Build( self, data );
- if ( self->temperatureField )
- Stg_Component_Build( self->temperatureField, data, False );
- if ( self->pressureField )
- Stg_Component_Build( self->pressureField, data, False );
+ if ( self->temperatureField )
+ Stg_Component_Build( self->temperatureField, data, False );
+ if ( self->pressureField )
+ Stg_Component_Build( self->pressureField, data, False );
- /* Sort out material extension stuff */
- self->materialExtHandle = Materials_Register_AddMaterialExtension(
- self->materials_Register,
- self->type,
- sizeof(BuoyancyForceTerm_MaterialExt) );
+ /* Sort out material extension stuff */
+ self->materialExtHandle=
+ Materials_Register_AddMaterialExtension(self->materials_Register,
+ self->type,
+ sizeof(BuoyancyForceTerm_MaterialExt));
- Journal_Firewall(densities.empty() && alphas.empty(),
- Journal_MyStream(Error_Type,self),
- "You may only create one instance of BuoyancyForceTerm");
+ Journal_Firewall(equations[0].empty() && equations[1].empty(),
+ Journal_MyStream(Error_Type,self),
+ "You may only create one instance of BuoyancyForceTerm");
- for ( material_I = 0 ; material_I < Materials_Register_GetCount( materials_Register ) ; material_I++) {
- material = Materials_Register_GetByIndex( materials_Register, material_I );
- materialExt = (BuoyancyForceTerm_MaterialExt*)ExtensionManager_GetFunc( material->extensionMgr, material, self->materialExtHandle );
+ for(material_I=0;
+ material_I<Materials_Register_GetCount(materials_Register);
+ material_I++)
+ {
+ material=Materials_Register_GetByIndex(materials_Register,material_I);
+ materialExt=(BuoyancyForceTerm_MaterialExt*)
+ ExtensionManager_GetFunc(material->extensionMgr,material,
+ self->materialExtHandle);
- /* Set the density */
- densities.push_back(Stg_ComponentFactory_GetString
- (cf,material->name,
- (Dictionary_Entry_Key)"densityEquation",""));
+ /* Set the density */
+ equations[0].push_back(Stg_ComponentFactory_GetString
+ (cf,material->name,
+ (Dictionary_Entry_Key)"densityEquation",
+ ""));
- materialExt->density =
- Stg_ComponentFactory_GetDouble
- (cf,material->name,(Dictionary_Entry_Key)"density",0.0);
+ materialExt->density =
+ Stg_ComponentFactory_GetDouble
+ (cf,material->name,(Dictionary_Entry_Key)"density",0.0);
- /* Set alpha */
- alphas.push_back(Stg_ComponentFactory_GetString
- (cf,material->name,
- (Dictionary_Entry_Key)"alphaEquation",""));
+ /* Set alpha */
+ equations[1].push_back(Stg_ComponentFactory_GetString
+ (cf,material->name,
+ (Dictionary_Entry_Key)"alphaEquation",
+ ""));
- materialExt->alpha =
- Stg_ComponentFactory_GetDouble
- (cf,material->name,(Dictionary_Entry_Key)"alpha",0.0);
- }
+ materialExt->alpha =
+ Stg_ComponentFactory_GetDouble
+ (cf,material->name,(Dictionary_Entry_Key)"alpha",0.0);
+ }
- /* Create Swarm Variables of each material swarm this ip swarm is mapped against */
- materialSwarms = IntegrationPointMapper_GetMaterialPointsSwarms( swarm->mapper, &(self->materialSwarmCount) );
- self->densitySwarmVariables = Memory_Alloc_Array( MaterialSwarmVariable*, self->materialSwarmCount, "DensityVariables" );
- self->alphaSwarmVariables = Memory_Alloc_Array( MaterialSwarmVariable*, self->materialSwarmCount, "AlphaVariables" );
+ /* Create Swarm Variables of each material swarm this ip swarm
+ is mapped against */
+ materialSwarms=IntegrationPointMapper_GetMaterialPointsSwarms
+ (swarm->mapper, &(self->materialSwarmCount) );
+ self->densitySwarmVariables=
+ Memory_Alloc_Array(MaterialSwarmVariable*, self->materialSwarmCount,
+ "DensityVariables" );
+ self->alphaSwarmVariables=
+ Memory_Alloc_Array(MaterialSwarmVariable*, self->materialSwarmCount,
+ "AlphaVariables");
- for ( materialSwarm_I = 0; materialSwarm_I < self->materialSwarmCount; ++materialSwarm_I ) {
- name = Stg_Object_AppendSuffix( materialSwarms[materialSwarm_I], (Name)"Density" );
- self->densitySwarmVariables[materialSwarm_I] = MaterialSwarmVariable_New(
- name,
- (AbstractContext*)self->context,
- materialSwarms[materialSwarm_I],
- 1,
- self->materials_Register,
- self->materialExtHandle,
- GetOffsetOfMember( *materialExt, density ) );
- Memory_Free( name );
+ for(materialSwarm_I=0; materialSwarm_I<self->materialSwarmCount;
+ ++materialSwarm_I )
+ {
+ name = Stg_Object_AppendSuffix( materialSwarms[materialSwarm_I], (Name)"Density" );
+ self->densitySwarmVariables[materialSwarm_I]=
+ MaterialSwarmVariable_New(name,
+ (AbstractContext*)self->context,
+ materialSwarms[materialSwarm_I],
+ 1,
+ self->materials_Register,
+ self->materialExtHandle,
+ GetOffsetOfMember( *materialExt, density ) );
+ Memory_Free( name );
- name = Stg_Object_AppendSuffix( materialSwarms[materialSwarm_I], (Name)"Alpha" );
- self->alphaSwarmVariables[materialSwarm_I] = MaterialSwarmVariable_New(
- name,
- (AbstractContext*)self->context,
- materialSwarms[materialSwarm_I],
- 1,
- self->materials_Register,
- self->materialExtHandle,
- GetOffsetOfMember( *materialExt, alpha ) );
- Memory_Free( name );
+ name=Stg_Object_AppendSuffix(materialSwarms[materialSwarm_I],
+ (Name)"Alpha");
+ self->alphaSwarmVariables[materialSwarm_I]=
+ MaterialSwarmVariable_New(name,
+ (AbstractContext*)self->context,
+ materialSwarms[materialSwarm_I],
+ 1,
+ self->materials_Register,
+ self->materialExtHandle,
+ GetOffsetOfMember( *materialExt, alpha ) );
+ Memory_Free( name );
- /* Build new Swarm Variables */
- Stg_Component_Build( self->densitySwarmVariables[materialSwarm_I], data, False );
- Stg_Component_Build( self->alphaSwarmVariables[materialSwarm_I], data, False );
- }
+ /* Build new Swarm Variables */
+ Stg_Component_Build(self->densitySwarmVariables[materialSwarm_I],
+ data, False);
+ Stg_Component_Build(self->alphaSwarmVariables[materialSwarm_I],
+ data, False);
+ }
}
-void _BuoyancyForceTerm_Initialise( void* forceTerm, void* data ) {
- BuoyancyForceTerm* self = (BuoyancyForceTerm*)forceTerm;
- Index i;
+void _BuoyancyForceTerm_Initialise( void* forceTerm, void* data )
+{
+ BuoyancyForceTerm* self=(BuoyancyForceTerm*)forceTerm;
+ Index i;
- _ForceTerm_Initialise( self, data );
+ _ForceTerm_Initialise( self, data );
- if ( self->temperatureField )
- Stg_Component_Initialise( self->temperatureField, data, False );
- if ( self->pressureField )
- Stg_Component_Initialise( self->pressureField, data, False );
+ if ( self->temperatureField )
+ Stg_Component_Initialise( self->temperatureField, data, False );
+ if ( self->pressureField )
+ Stg_Component_Initialise( self->pressureField, data, False );
- for ( i = 0; i < self->materialSwarmCount; ++i ) {
- Stg_Component_Initialise( self->densitySwarmVariables[i], data, False );
- Stg_Component_Initialise( self->alphaSwarmVariables[i], data, False );
- }
+ for ( i = 0; i < self->materialSwarmCount; ++i ) {
+ Stg_Component_Initialise( self->densitySwarmVariables[i], data, False );
+ Stg_Component_Initialise( self->alphaSwarmVariables[i], data, False );
+ }
}
-void _BuoyancyForceTerm_Execute( void* forceTerm, void* data ) {
- BuoyancyForceTerm* self = (BuoyancyForceTerm*)forceTerm;
-
- _ForceTerm_Execute( self, data );
+void _BuoyancyForceTerm_Execute( void* forceTerm, void* data )
+{
+ BuoyancyForceTerm* self = (BuoyancyForceTerm*)forceTerm;
+ _ForceTerm_Execute( self, data );
}
-void _BuoyancyForceTerm_Destroy( void* forceTerm, void* data ) {
- BuoyancyForceTerm* self = (BuoyancyForceTerm*)forceTerm;
- Index i;
+void _BuoyancyForceTerm_Destroy( void* forceTerm, void* data )
+{
+ BuoyancyForceTerm* self = (BuoyancyForceTerm*)forceTerm;
+ Index i;
- for ( i = 0; i < self->materialSwarmCount; ++i ) {
- _Stg_Component_Delete( self->densitySwarmVariables[i] );
- _Stg_Component_Delete( self->alphaSwarmVariables[i] );
- }
+ for ( i = 0; i < self->materialSwarmCount; ++i ) {
+ _Stg_Component_Delete( self->densitySwarmVariables[i] );
+ _Stg_Component_Delete( self->alphaSwarmVariables[i] );
+ }
- FreeArray( self->gHat );
+ FreeArray( self->gHat );
- Memory_Free( self->densitySwarmVariables );
- Memory_Free( self->alphaSwarmVariables );
+ Memory_Free( self->densitySwarmVariables );
+ Memory_Free( self->alphaSwarmVariables );
- _ForceTerm_Destroy( forceTerm, data );
+ _ForceTerm_Destroy( forceTerm, data );
}
void _BuoyancyForceTerm_AssembleElement(void* forceTerm,
@@ -365,43 +417,16 @@ void _BuoyancyForceTerm_AssembleElement(
double* elForceVec)
{
BuoyancyForceTerm* self = (BuoyancyForceTerm*) forceTerm;
- IntegrationPoint* particle;
- Particle_InCellIndex cParticle_I;
- Particle_InCellIndex cellParticleCount;
- Element_NodeIndex elementNodeCount;
Dimension_Index dim = forceVector->dim;
IntegrationPointsSwarm* swarm=(IntegrationPointsSwarm*)self->integrationSwarm;
FeMesh* mesh = forceVector->feVariable->feMesh;
- Node_ElementLocalIndex eNode_I;
- Cell_Index cell_I;
- ElementType* elementType;
- Dof_Index nodeDofCount;
- double gravity;
- double detJac = 0.0;
- double factor;
- double Ni[27];
- double force;
- double* xi;
- FeVariable* temperatureField = self->temperatureField;
- double temperature = 0.0;
- FeVariable* pressureField = self->pressureField;
- double pressure = 0.0;
- double background_density = 0.0;
- double* gHat;
- unsigned d_i;
- double totalWeight = 0.0;
- double adjustFactor = 0.0;
- double density(0);
- double alpha(0);
- std::string density_equation, alpha_equation;
-
- elementType = FeMesh_GetElementType( mesh, lElement_I );
- elementNodeCount = elementType->nodeCount;
- nodeDofCount = dim;
- cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
- cellParticleCount = swarm->cellParticleCountTbl[cell_I];
- gHat = self->gHat;
+ ElementType* elementType = FeMesh_GetElementType( mesh, lElement_I );
+ Element_NodeIndex elementNodeCount = elementType->nodeCount;
+ Dof_Index nodeDofCount = dim;
+ Cell_Index cell_I=
+ CellLayout_MapElementIdToCellId(swarm->cellLayout,lElement_I);
+ Particle_InCellIndex cellParticleCount = swarm->cellParticleCountTbl[cell_I];
/* adjust & adjustFactor -- 20060411 Alan
*
@@ -413,12 +438,15 @@ void _BuoyancyForceTerm_AssembleElement(
* not cover the whole domain (ie - I used it to do dave.m's test of
* 1 swarm for blob, 1 swarm for background)
*/
- if (self->adjust)
+ double adjustFactor;
+ if(self->adjust)
{
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ )
+ double totalWeight = 0.0;
+ for(Particle_InCellIndex cParticle_I=0; cParticle_I<cellParticleCount;
+ cParticle_I++)
{
- particle=(IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,
- cParticle_I);
+ IntegrationPoint* particle=
+ (IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,cParticle_I);
totalWeight += particle->weight;
}
adjustFactor = swarm->weights->cellLocalVolume / totalWeight;
@@ -429,153 +457,179 @@ void _BuoyancyForceTerm_AssembleElement(
}
/* Use an average density and alpha for OneToMany mapper */
+ double density(0), alpha(0);
bool one_to_many=Stg_Class_IsInstance(swarm->mapper,OneToManyMapper_Type);
if(one_to_many)
{
- OneToManyMapper *mapper=(OneToManyMapper*)(swarm->mapper);
- IntegrationPointsSwarm* OneToMany_swarm=mapper->swarm;
- int OneToMany_cell=
- CellLayout_MapElementIdToCellId(OneToMany_swarm->cellLayout,
- lElement_I);
- int num_particles=OneToMany_swarm->cellParticleCountTbl[OneToMany_cell];
-
- double weight(0);
- for(int ii=0;ii<num_particles;++ii)
- {
- IntegrationPoint *OneToMany_particle=
- (IntegrationPoint*)Swarm_ParticleInCellAt(OneToMany_swarm,
- OneToMany_cell,ii);
- weight+=OneToMany_particle->weight;
- double ttmp=IntegrationPointMapper_GetDoubleFromMaterial
- (OneToMany_swarm->mapper, OneToMany_particle,
- self->materialExtHandle,
- offsetof(BuoyancyForceTerm_MaterialExt, density));
- density+=ttmp*OneToMany_particle->weight;
-
- ttmp=IntegrationPointMapper_GetDoubleFromMaterial
- (OneToMany_swarm->mapper, OneToMany_particle,
- self->materialExtHandle,
- offsetof(BuoyancyForceTerm_MaterialExt, alpha));
- alpha+=ttmp*OneToMany_particle->weight;
- }
- density/=weight;
- alpha/=weight;
+ BuoyancyForceTerm_average_density_alpha(self,swarm,mesh,lElement_I,
+ &density,&alpha);
}
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- Coord coord;
- particle=(IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,
- cParticle_I);
- xi=particle->xi;
+ for(Particle_InCellIndex cParticle_I=0; cParticle_I<cellParticleCount;
+ cParticle_I++)
+ {
+ IntegrationPoint* particle=
+ (IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,cParticle_I);
+ double *xi=particle->xi;
- detJac=ElementType_JacobianDeterminant(elementType,mesh,lElement_I,xi,dim);
- ElementType_EvaluateShapeFunctionsAt(elementType,xi,Ni);
+ double detJac=ElementType_JacobianDeterminant(elementType,mesh,
+ lElement_I,xi,dim);
+ double Ni[27];
+ ElementType_EvaluateShapeFunctionsAt(elementType,xi,Ni);
- /* Get parameters */
- if(temperatureField)
- FeVariable_InterpolateFromMeshLocalCoord
- (temperatureField, mesh, lElement_I, xi, &temperature);
+ /* Get parameters */
+ double temperature(0), pressure(0);
+ if(self->temperatureField)
+ FeVariable_InterpolateFromMeshLocalCoord
+ (self->temperatureField, mesh, lElement_I, xi, &temperature);
- if(pressureField)
- FeVariable_InterpolateFromMeshLocalCoord
- (pressureField, mesh, lElement_I, xi, &pressure);
+ if(self->pressureField)
+ FeVariable_InterpolateFromMeshLocalCoord
+ (self->pressureField, mesh, lElement_I, xi, &pressure);
- if(self->hydrostaticTerm)
- {
- FeMesh_CoordLocalToGlobal(mesh, cell_I, xi, coord);
- background_density=
- HydrostaticTerm_Density(self->hydrostaticTerm,coord);
- if(pressureField)
- pressure+=HydrostaticTerm_Pressure(self->hydrostaticTerm,
- coord);
- }
+ double background_density = 0.0;
+ if(self->hydrostaticTerm)
+ {
+ Coord coord;
+ FeMesh_CoordLocalToGlobal(mesh, cell_I, xi, coord);
+ background_density=
+ HydrostaticTerm_Density(self->hydrostaticTerm,coord);
+ if(self->pressureField)
+ pressure+=HydrostaticTerm_Pressure(self->hydrostaticTerm,
+ coord);
+ }
- if(!one_to_many)
- {
- /* Handle case where we are using gauss swarms with
- NearestNeighborMapper instead of a material
- swarm */
- IntegrationPointsSwarm* NNswarm(swarm);
- IntegrationPoint* NNparticle(particle);
- NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
+ if(!one_to_many)
+ {
+ /* Handle case where we are using gauss swarms with
+ NearestNeighborMapper instead of a material
+ swarm */
+ IntegrationPointsSwarm* NNswarm(swarm);
+ IntegrationPoint* NNparticle(particle);
+ NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
- /* Get density and alpha */
- unsigned int material_index=
- IntegrationPointMapper_GetMaterialIndexOn(NNswarm->mapper,
- NNparticle);
- Journal_Firewall(material_index<densities.size()
- && material_index<alphas.size(),
- Journal_MyStream(Error_Type,self),
- "Bad material index %d",material_index);
+ /* Get density and alpha */
+ density=BuoyancyForceTerm_density(self,mesh,lElement_I,cell_I,
+ NNswarm->mapper,NNparticle);
+ alpha=BuoyancyForceTerm_alpha(self,mesh,lElement_I,cell_I,
+ NNswarm->mapper,NNparticle);
+ }
+ /* Calculate Force */
+ double gravity[dim];
+ BuoyancyForceTerm_CalcGravity(self,mesh,lElement_I,xi,gravity);
- density_equation = densities[material_index];
- if(!density_equation.empty())
- {
- std::stringstream ss;
+ double rho=(density*(1.0-alpha*temperature) - background_density);
+ double factor = detJac * particle->weight * adjustFactor * rho;
- ss << "T=" << temperature
- << "; p=" << pressure
- << "; " << density_equation;
- density=Equation_eval(coord,
- (DomainContext *)(self->context),
- ss.str());
- }
- else
- {
- density = IntegrationPointMapper_GetDoubleFromMaterial
- (NNswarm->mapper, NNparticle, self->materialExtHandle,
- offsetof(BuoyancyForceTerm_MaterialExt, density));
- }
-
- alpha_equation = alphas[material_index];
- if(!alpha_equation.empty())
- {
- std::stringstream ss;
-
- ss << "T=" << temperature
- << "; p=" << pressure
- << "; " << density_equation;
- alpha=Equation_eval(coord,
- (DomainContext *)(self->context),
- ss.str());
- }
- else
- {
- alpha = IntegrationPointMapper_GetDoubleFromMaterial
- (NNswarm->mapper, NNparticle, self->materialExtHandle,
- offsetof(BuoyancyForceTerm_MaterialExt, alpha));
- }
- }
- /* Calculate Force */
- gravity = BuoyancyForceTerm_CalcGravity(self,(Swarm*)swarm,lElement_I,
- particle);
- force=gravity*(density*(1.0-alpha*temperature)
- - background_density);
- factor = detJac * particle->weight * adjustFactor * force;
-
- /* Apply force in the correct direction */
- for(eNode_I=0; eNode_I<elementNodeCount; eNode_I++)
- {
- if(gHat)
- {
- for(d_i=0; d_i<dim; d_i++)
- elForceVec[ eNode_I * nodeDofCount + d_i ] +=
- gHat[d_i] * factor * Ni[ eNode_I ] ;
- }
- else
- {
- elForceVec[ eNode_I * nodeDofCount + J_AXIS ] +=
- -1.0 * factor * Ni[ eNode_I ] ;
- }
- }
- }
+ /* Apply force in the correct direction */
+ for(Node_ElementLocalIndex eNode_I=0; eNode_I<elementNodeCount; eNode_I++)
+ {
+ for(unsigned d_i=0; d_i<dim; d_i++)
+ elForceVec[ eNode_I * nodeDofCount + d_i ] +=
+ gravity[d_i] * factor * Ni[ eNode_I ] ;
+ }
+ }
}
/* The default implementation is for the gravity to be constant. */
-double _BuoyancyForceTerm_CalcGravity( void* forceTerm, Swarm* swarm, Element_DomainIndex dElement_I, void* particle ) {
- BuoyancyForceTerm* self = (BuoyancyForceTerm*) forceTerm;
-
- return self->gravity;
+void _BuoyancyForceTerm_CalcGravity(BuoyancyForceTerm *self, FeMesh *mesh,
+ Element_DomainIndex dElement_I,
+ double* xi, double *gravity)
+{
+ for(unsigned int d=0;d<mesh->generator->nDims;++d)
+ gravity[d]=self->gHat[d];
}
+void BuoyancyForceTerm_average_density_alpha(BuoyancyForceTerm *self,
+ IntegrationPointsSwarm* input_swarm,
+ FeMesh* mesh,
+ const Element_LocalIndex &lElement_I,
+ double *density, double *alpha)
+{
+ *density=*alpha=0;
+ IntegrationPointsSwarm* swarm(input_swarm);
+ bool one_to_many=Stg_Class_IsInstance(swarm->mapper,OneToManyMapper_Type);
+ if(one_to_many)
+ {
+ swarm=((OneToManyMapper*)(swarm->mapper))->swarm;
+ }
+ int cell_I=CellLayout_MapElementIdToCellId(swarm->cellLayout,
+ lElement_I);
+ int num_particles=swarm->cellParticleCountTbl[cell_I];
+ double weight(0);
+ for(int ii=0;ii<num_particles;++ii)
+ {
+ IntegrationPoint *particle=
+ (IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,ii);
+
+ weight+=particle->weight;
+ *density+=
+ BuoyancyForceTerm_density(self,mesh,lElement_I,cell_I,
+ swarm->mapper,particle)*particle->weight;
+
+ *alpha+=BuoyancyForceTerm_alpha(self,mesh,lElement_I,cell_I,
+ swarm->mapper,particle)*particle->weight;
+
+ }
+ *density/=weight;
+ *alpha/=weight;
+}
+
+double BuoyancyForceTerm_density_alpha(BuoyancyForceTerm *self,
+ FeMesh* mesh,
+ const Element_LocalIndex &lElement_I,
+ const Cell_Index &cell_I,
+ IntegrationPointMapper* mapper,
+ IntegrationPoint* particle,
+ const int &type)
+{
+ unsigned int material_index=
+ IntegrationPointMapper_GetMaterialIndexOn(mapper,particle);
+
+ Journal_Firewall(material_index<equations[type].size(),
+ Journal_MyStream(Error_Type,self),
+ "Bad material index %d",material_index);
+
+ std::string equation = equations[type][material_index];
+ double result;
+ if(!equation.empty())
+ {
+ std::stringstream ss;
+
+ double temperature(0), pressure(0);
+ if(self->temperatureField)
+ FeVariable_InterpolateFromMeshLocalCoord
+ (self->temperatureField, mesh, lElement_I, particle->xi, &temperature);
+
+ if(self->pressureField)
+ FeVariable_InterpolateFromMeshLocalCoord
+ (self->pressureField, mesh, lElement_I, particle->xi, &pressure);
+
+ Coord coord;
+ FeMesh_CoordLocalToGlobal(mesh, cell_I, particle->xi, coord);
+
+ ss << "T=" << temperature
+ << "; p=" << pressure
+ << "; " << equation;
+ result=Equation_eval(coord,
+ (DomainContext *)(self->context),
+ ss.str());
+ }
+ else
+ {
+ if(type==0)
+ {
+ result=IntegrationPointMapper_GetDoubleFromMaterial
+ (mapper, particle, self->materialExtHandle,
+ offsetof(BuoyancyForceTerm_MaterialExt, density));
+ }
+ else
+ {
+ result=IntegrationPointMapper_GetDoubleFromMaterial
+ (mapper, particle, self->materialExtHandle,
+ offsetof(BuoyancyForceTerm_MaterialExt, alpha));
+ }
+ }
+ return result;
+}
diff -r 30411e58fae6 -r c1e2c18c39f8 Utils/src/BuoyancyForceTerm.h
--- a/Utils/src/BuoyancyForceTerm.h Sat Apr 07 15:10:18 2012 -0700
+++ b/Utils/src/BuoyancyForceTerm.h Mon Apr 30 13:50:43 2012 -0700
@@ -44,7 +44,7 @@
#ifndef __PICellerator_Utils_BuoyancyForceTerm_h__
#define __PICellerator_Utils_BuoyancyForceTerm_h__
- typedef double (BuoyancyForceTerm_CalcGravityFunction) ( void* forceTerm, Swarm* swarm, Element_LocalIndex lElement_I, void* particle );
+ typedef void (BuoyancyForceTerm_CalcGravityFunction) (BuoyancyForceTerm *self, FeMesh *mesh, Element_DomainIndex dElement_I, double* xi, double *gravity);
/** Textual name of this class */
extern const Type BuoyancyForceTerm_Type;
@@ -70,6 +70,7 @@
double gravity; \
double* gHat; \
Bool adjust; \
+ Bool damping; \
Materials_Register* materials_Register; \
ExtensionInfo_Index materialExtHandle; \
MaterialSwarmVariable** densitySwarmVariables; \
@@ -89,7 +90,9 @@
FeVariable* temperatureField,
FeVariable* pressureField,
double gravity,
+ double* gHat,
Bool adjust,
+ Bool damping,
Materials_Register* materials_Register,
HydrostaticTerm* hydrostaticTerm);
@@ -113,7 +116,9 @@
FeVariable* temperatureField,
FeVariable* pressureField,
double gravity,
+ double* gHat,
Bool adjust,
+ Bool damping,
Materials_Register* materials_Register,
HydrostaticTerm* hydrostaticTerm );
@@ -135,10 +140,49 @@
void _BuoyancyForceTerm_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) ;
- double _BuoyancyForceTerm_CalcGravity( void* forceTerm, Swarm* swarm, Element_LocalIndex lElement_I, void* particle ) ;
+ void _BuoyancyForceTerm_CalcGravity(BuoyancyForceTerm *self,
+ FeMesh *mesh,
+ Element_DomainIndex dElement_I,
+ double* xi, double *gravity);
- #define BuoyancyForceTerm_CalcGravity( forceTerm, swarm, lElement_I, particle ) \
- (( (BuoyancyForceTerm*) forceTerm )->_calcGravity( forceTerm, swarm, lElement_I, particle ) )
+ #define BuoyancyForceTerm_CalcGravity( forceTerm, mesh, dElement_I, xi, gravity ) \
+ (( (BuoyancyForceTerm*) forceTerm )->_calcGravity( forceTerm, mesh, dElement_I, xi, gravity ) )
+
+ void BuoyancyForceTerm_average_density_alpha(BuoyancyForceTerm *self,
+ IntegrationPointsSwarm* input_swarm,
+ FeMesh* mesh,
+ const Element_LocalIndex &lElement_I,
+ double *density, double *alpha);
+
+ double BuoyancyForceTerm_density_alpha(BuoyancyForceTerm *self,
+ FeMesh* mesh,
+ const Element_LocalIndex &lElement_I,
+ const Cell_Index &cell_I,
+ IntegrationPointMapper* mapper,
+ IntegrationPoint* particle,
+ const int &type);
+
+ inline double BuoyancyForceTerm_density(BuoyancyForceTerm *self,
+ FeMesh* mesh,
+ const Element_LocalIndex &lElement_I,
+ const Cell_Index &cell_I,
+ IntegrationPointMapper* mapper,
+ IntegrationPoint* particle)
+ {
+ return BuoyancyForceTerm_density_alpha(self,mesh,lElement_I,cell_I,
+ mapper,particle,0);
+ }
+
+ inline double BuoyancyForceTerm_alpha(BuoyancyForceTerm *self,
+ FeMesh* mesh,
+ const Element_LocalIndex &lElement_I,
+ const Cell_Index &cell_I,
+ IntegrationPointMapper* mapper,
+ IntegrationPoint* particle)
+ {
+ return BuoyancyForceTerm_density_alpha(self,mesh,lElement_I,cell_I,
+ mapper,particle,1);
+ }
#endif
diff -r 30411e58fae6 -r c1e2c18c39f8 Utils/src/Init.cxx
--- a/Utils/src/Init.cxx Sat Apr 07 15:10:18 2012 -0700
+++ b/Utils/src/Init.cxx Mon Apr 30 13:50:43 2012 -0700
@@ -62,11 +62,13 @@ Bool PICellerator_Utils_Init( int* argc,
Stg_ComponentRegister_Add( componentsRegister, BuoyancyForceTerm_Type, (Name)"0", _BuoyancyForceTerm_DefaultNew );
Stg_ComponentRegister_Add( componentsRegister, BuoyancyForceTermThermoChem_Type, (Name)"0", _BuoyancyForceTermThermoChem_DefaultNew );
Stg_ComponentRegister_Add( componentsRegister, DiffusionSMT_Type, (Name)"0", _DiffusionSMT_DefaultNew );
+ Stg_ComponentRegister_Add( componentsRegister, BuoyancyDampingTerm_Type, (Name)"0", _BuoyancyDampingTerm_DefaultNew );
Stg_ComponentRegister_Add( componentsRegister, HydrostaticTerm_Type, (Name)"0", _HydrostaticTerm_DefaultNew );
RegisterParent( BuoyancyForceTerm_Type, ForceTerm_Type );
RegisterParent( BuoyancyForceTermThermoChem_Type, ForceTerm_Type );
RegisterParent( DiffusionSMT_Type, StiffnessMatrixTerm_Type );
+ RegisterParent( BuoyancyDampingTerm_Type, StiffnessMatrixTerm_Type );
RegisterParent( MaterialSwarmVariable_Type, SwarmVariable_Type );
RegisterParent( PCDVC_Type, DVCWeights_Type );
RegisterParent( HydrostaticTerm_Type, Stg_Component_Type );
diff -r 30411e58fae6 -r c1e2c18c39f8 Utils/src/Utils.h
--- a/Utils/src/Utils.h Sat Apr 07 15:10:18 2012 -0700
+++ b/Utils/src/Utils.h Mon Apr 30 13:50:43 2012 -0700
@@ -46,6 +46,7 @@
#include "types.h"
#include "BuoyancyForceTerm.h"
+ #include "BuoyancyDampingTerm.h"
#include "BuoyancyForceTermThermoChem.h"
#include "DiffusionSMT.h"
#include "MaterialSwarmVariable.h"
diff -r 30411e58fae6 -r c1e2c18c39f8 Utils/src/types.h
--- a/Utils/src/types.h Sat Apr 07 15:10:18 2012 -0700
+++ b/Utils/src/types.h Mon Apr 30 13:50:43 2012 -0700
@@ -54,9 +54,10 @@
#ifndef __PICellerator_Utils_types_h__
#define __PICellerator_Utils_types_h__
- typedef struct MaterialSwarmVariable MaterialSwarmVariable;
- typedef struct BuoyancyForceTerm BuoyancyForceTerm;
- typedef struct BuoyancyForceTermThermoChem BuoyancyForceTermThermoChem;
+ typedef struct MaterialSwarmVariable MaterialSwarmVariable;
+ typedef struct BuoyancyForceTerm BuoyancyForceTerm;
+ typedef struct BuoyancyDampingTerm BuoyancyDampingTerm;
+ typedef struct BuoyancyForceTermThermoChem BuoyancyForceTermThermoChem;
typedef struct DiffusionSMT DiffusionSMT;
typedef struct HydrostaticTerm HydrostaticTerm;
typedef struct PCDVC PCDVC;
More information about the CIG-COMMITS
mailing list