[cig-commits] r5677 - in long/3D/Gale/trunk/src/StgFEM: . Discretisation/src Discretisation/tests

walter at geodynamics.org walter at geodynamics.org
Fri Jan 5 11:37:32 PST 2007


Author: walter
Date: 2007-01-05 11:37:31 -0800 (Fri, 05 Jan 2007)
New Revision: 5677

Modified:
   long/3D/Gale/trunk/src/StgFEM/
   long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeEquationNumber.c
   long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c
   long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testFeEquationNumber-LinkedDofs.c
Log:
 r909 at earth (orig r688):  LukeHodkinson | 2007-01-02 21:31:16 -0800
 Added linked degrees of freedom back in.
 
 



Property changes on: long/3D/Gale/trunk/src/StgFEM
___________________________________________________________________
Name: svk:merge
   - 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:687
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
   + 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:688
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669

Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeEquationNumber.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeEquationNumber.c	2007-01-05 19:37:28 UTC (rev 5676)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeEquationNumber.c	2007-01-05 19:37:31 UTC (rev 5677)
@@ -1971,10 +1971,9 @@
 	Node_GlobalIndex gNode_I;
 	Node_GlobalIndex nodeGlobalCount = FeMesh_GetNodeGlobalSize( feMesh );
 
+	MPI_Comm_rank( comm, (int*)&rank );
 	Journal_Printf( stream, "%d: *** Printing destination array ***\n", rank );
 
-	MPI_Comm_rank( comm, (int*)&rank );
-
 	for (gNode_I =0; gNode_I < nodeGlobalCount; gNode_I++) {
 		Node_DomainIndex dNode_I;
 
@@ -2026,7 +2025,11 @@
 			/* print the nodes and dofs */
 			for ( elLocalNode_I = 0; elLocalNode_I < numNodesAtElement; elLocalNode_I++ ) {
 				/* look up processor local node number. */
-				Element_LocalIndex currNode = incNodes[elLocalNode_I];
+				Element_LocalIndex currNode = incNodes[elLocalNode_I == 2 ? 3 : 
+								       elLocalNode_I == 3 ? 2 : 
+								       elLocalNode_I == 6 ? 7 : 
+								       elLocalNode_I == 7 ? 6 : 
+								       elLocalNode_I];
 				/* get the number of dofs at current node */
 				Dof_Index currNodeNumDofs = self->dofLayout->dofCounts[ currNode ];
 				Dof_Index nodeLocalDof_I;
@@ -2267,6 +2270,7 @@
 	unsigned		maxDofs;
 	unsigned*		tuples;
 	Decomp_Sync_Array*	array;
+	LinkedDofInfo*		links;
 	unsigned		e_i, n_i, dof_i;
 
 	assert( self );
@@ -2282,11 +2286,21 @@
 	nDomainEls = FeMesh_GetElementDomainSize( feMesh );
 	nLocalNodes = FeMesh_GetNodeLocalSize( feMesh );
 	nNodalDofs = self->dofLayout->dofCounts;
+	links = self->linkedDofInfo;
 
 	/* Allocate for destination array. */
 	dstArray = Memory_Alloc_2DComplex( int, nDomainNodes, nNodalDofs, 
 					   "FeEquationNumber::destinationArray" );
 
+	/* If needed, allocate for linked equation numbers. */
+	if( links ) {
+		unsigned	s_i;
+
+		links->eqNumsOfLinkedDofs = ReallocArray( links->eqNumsOfLinkedDofs, int, links->linkedDofSetsCount );
+		for( s_i = 0; s_i < links->linkedDofSetsCount; s_i++ )
+			links->eqNumsOfLinkedDofs[s_i] = -1;
+	}
+
 	/* Allocate for the location matrix. */
 	nLocMatDofs = NULL;
 	locMat = AllocArray( int**, nDomainEls );
@@ -2308,8 +2322,15 @@
 
 		for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
 			varInd = self->dofLayout->varIndices[n_i][dof_i];
-			if( !self->bcs || !VariableCondition_IsCondition( self->bcs, n_i, varInd ) )
-				dstArray[n_i][dof_i] = curEqNum++;
+			if( !self->bcs || !VariableCondition_IsCondition( self->bcs, n_i, varInd ) ) {
+				if( links && links->linkedDofTbl[n_i][dof_i] != -1 ) {
+					if( links->eqNumsOfLinkedDofs[links->linkedDofTbl[n_i][dof_i]] == -1 )
+						links->eqNumsOfLinkedDofs[links->linkedDofTbl[n_i][dof_i]] = curEqNum++;
+					dstArray[n_i][dof_i] = links->eqNumsOfLinkedDofs[links->linkedDofTbl[n_i][dof_i]];
+				}
+				else
+					dstArray[n_i][dof_i] = curEqNum++;
+			}
 			else
 				dstArray[n_i][dof_i] = -1;
 		}
@@ -2325,13 +2346,31 @@
 	if( rank < nProcs - 1 )
 		MPI_Send( &subTotal, 1, MPI_UNSIGNED, rank + 1, 6669, comm );
 
+	/* Reduce to find lowest linked DOFs. */
+	if( links ) {
+		unsigned	lowest;
+		unsigned	s_i;
+
+		for( s_i = 0; s_i < links->linkedDofSetsCount; s_i++ ) {
+			if( links->eqNumsOfLinkedDofs[s_i] != -1 )
+				links->eqNumsOfLinkedDofs[s_i] += base;
+			MPI_Allreduce( links->eqNumsOfLinkedDofs + s_i, &lowest, 1, MPI_UNSIGNED, MPI_MIN, comm );
+			links->eqNumsOfLinkedDofs[s_i] = lowest;
+			assert( links->eqNumsOfLinkedDofs[s_i] >= 0 );
+		}
+	}
+
 	/* Modify existing destination array and dump to a tuple array. */
 	tuples = AllocArray( unsigned, nDomainNodes * maxDofs );
 	for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
 		for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
 			varInd = self->dofLayout->varIndices[n_i][dof_i];
-			if( !self->bcs || !VariableCondition_IsCondition( self->bcs, n_i, varInd ) )
-				dstArray[n_i][dof_i] += base;
+			if( !self->bcs || !VariableCondition_IsCondition( self->bcs, n_i, varInd ) ) {
+				if( links && links->linkedDofTbl[n_i][dof_i] != -1 )
+					dstArray[n_i][dof_i] = links->eqNumsOfLinkedDofs[links->linkedDofTbl[n_i][dof_i]];
+				else
+					dstArray[n_i][dof_i] += base;
+			}
 			tuples[n_i * maxDofs + dof_i] = dstArray[n_i][dof_i];
 		}
 	}

Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c	2007-01-05 19:37:28 UTC (rev 5676)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c	2007-01-05 19:37:31 UTC (rev 5677)
@@ -811,11 +811,11 @@
 {
 	FeVariable*		self = (FeVariable*)feVariable;
 	InterpolationResult	retValue;
-	unsigned		elDim, elInd;
+	unsigned		elInd;
 
 	/* locate which mesh element given coord is in : use inclusive upper boundaries to save
 	   the need to use shadow space if possible */
-	if( !Mesh_Search( self->feMesh, globalCoord, &elDim, &elInd ) ) {
+	if( !Mesh_SearchElements( self->feMesh, globalCoord, &elInd ) ) {
 		Bool			outsideGlobal = False;
 		double			min[3], max[3];
 		Dimension_Index		dim_I=0;
@@ -836,17 +836,8 @@
 	}
 	else /* We found the coord is within a local or shadow element */ {
 		ElementType*		elementType = NULL;
-		unsigned		nInc, *inc;
 
-		if( elDim != Mesh_GetDimSize( self->feMesh ) ) {
-			Mesh_GetIncidence( self->feMesh, elDim, elInd, Mesh_GetDimSize( self->feMesh ), 
-					   &nInc, &inc );
-			assert( nInc );
-			*elementCoordInPtr = inc[0];
-		}
-		else
-			*elementCoordInPtr = elInd;
-	
+		*elementCoordInPtr = elInd;
 		if ( elInd < FeMesh_GetElementLocalSize( self->feMesh ) ) {
 			retValue = LOCAL;
 		}
@@ -1081,7 +1072,6 @@
 
 void FeVariable_InterpolateDerivatives_WithGNx( void* _feVariable, Element_LocalIndex lElement_I, double** GNx, double* value ) {
 	FeVariable*             self        = (FeVariable*) _feVariable;
-	ElementType*            elementType = FeMesh_GetElementType( self->feMesh, lElement_I );
 	Node_ElementLocalIndex  elLocalNode_I;
 	Node_LocalIndex         lNode_I;
 	Dof_Index               dof_I;
@@ -1491,9 +1481,8 @@
 	FeVariable*          self               = (FeVariable*)         feVariable;
 	Swarm*               swarm              = (Swarm*)              _swarm;
 	FeMesh*			feMesh             = self->feMesh;
-	FeMesh*			geometryMesh;
+	FeMesh*			mesh;
 	ElementType*         elementType;
-	ElementType*         geometryElementType;
 	Cell_LocalIndex      cell_I;
 	Particle_InCellIndex cParticle_I;
 	Particle_InCellIndex cellParticleCount;
@@ -1506,9 +1495,8 @@
 	integral = 0.0;
 
 	/* Use feVariable's mesh as geometry mesh if one isn't passed in */
-	geometryMesh = (feMesh->elMesh == (Mesh*)feMesh) ? feMesh : feMesh->elMesh;
-	elementType = FeMesh_GetElementType( geometryMesh, dElement_I );
-	geometryElementType = FeMesh_GetElementType( geometryMesh, dElement_I );
+	mesh = (feMesh->elMesh ? (FeMesh*)feMesh->elMesh : feMesh);
+	elementType = FeMesh_GetElementType( mesh, dElement_I );
 
 	/* Determine number of particles in element */
 	cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, dElement_I );
@@ -1524,7 +1512,7 @@
 
 		/* Calculate Determinant of Jacobian */
 		detJac = ElementType_JacobianDeterminant_AxisIndependent( 
-				geometryElementType, geometryMesh, dElement_I, particle->xi, dim, axis0, axis1, axis2 );
+				elementType, mesh, dElement_I, particle->xi, dim, axis0, axis1, axis2 );
 
 		/* Sum Integral */
 		integral += detJac * particle->weight * value;
@@ -1577,9 +1565,11 @@
 	Dimension_Index            dim                = self->dim;
 	double                     integral;
 	double                     layerThickness     = 0.0;
+	double			sendThickness;
 	Grid*			vertGrid;
 	unsigned*		inds;
 	double			heights[2];
+	unsigned		localInd[2], globalInd[2];
 	double			*min, *max;
 	unsigned		d_i;
 	
@@ -1595,16 +1585,26 @@
 		else
 			inds[d_i] = layerIndex;
 	}
-	heights[0] = Mesh_GetVertex( self->feMesh, Grid_Project( vertGrid, inds ) )[layerAxis];
+	globalInd[0] = Grid_Project( vertGrid, inds );
 	inds[layerAxis]++;
-	heights[1] = Mesh_GetVertex( self->feMesh, Grid_Project( vertGrid, inds ) )[layerAxis];
+	globalInd[1] = Grid_Project( vertGrid, inds );
+	if( Mesh_GlobalToDomain( self->feMesh, MT_VERTEX, globalInd[0], &localInd[0] ) && 
+	    Mesh_GlobalToDomain( self->feMesh, MT_VERTEX, globalInd[1], &localInd[1] ) )
+	{
+		heights[0] = Mesh_GetVertex( self->feMesh, localInd[0] )[layerAxis];
+		heights[1] = Mesh_GetVertex( self->feMesh, localInd[1] )[layerAxis];
+		sendThickness = heights[1] - heights[0];
+	}
+	else {
+		sendThickness = 0.0;
+	}
+	MPI_Allreduce( &sendThickness, &layerThickness, 1, MPI_DOUBLE, MPI_MAX, self->communicator );
 	FreeArray( inds );
-	layerThickness = heights[1] - heights[0];
 
 	min = Memory_Alloc_Array_Unnamed( double, Mesh_GetDimSize( self->feMesh ) );
 	max = Memory_Alloc_Array_Unnamed( double, Mesh_GetDimSize( self->feMesh ) );
 	Mesh_GetGlobalCoordRange( self->feMesh, min, max );
-	integral /= layerThickness * ( max[ aAxis ] - min[ aAxis ] );
+	integral /= layerThickness * (max[aAxis] - min[aAxis]);
 	if ( dim == 3 )
 		integral /= max[ bAxis ] - min[ bAxis ];
 	FreeArray( min );
@@ -1642,9 +1642,11 @@
 			continue;
 
 		/* Check if element is local */
-		insist( FeMesh_ElementGlobalToDomain( self->feMesh, gElement_I, &lElement_I ) );
-		if( lElement_I >= FeMesh_GetElementLocalSize( self->feMesh ) )
+		if( !FeMesh_ElementGlobalToDomain( self->feMesh, gElement_I, &lElement_I ) || 
+		    lElement_I >= FeMesh_GetElementLocalSize( self->feMesh ) )
+		{
 			continue;
+		}
 
 		elementIntegral = FeVariable_IntegrateElement_AxisIndependent( self, swarm, lElement_I, dim, axis0, axis1, axis2 );
 		Journal_DPrintfL( self->debug, 2, "Integral of element %d was %f\n", lElement_I, elementIntegral );

Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testFeEquationNumber-LinkedDofs.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testFeEquationNumber-LinkedDofs.c	2007-01-05 19:37:28 UTC (rev 5676)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testFeEquationNumber-LinkedDofs.c	2007-01-05 19:37:31 UTC (rev 5677)
@@ -63,10 +63,10 @@
 void _SetDt( void* context, double dt ) {
 }
 
-void Test_FeEquationNumberRun_Regular( Dictionary* dictionary, void* context, const char* nLayoutType, IJK elSizes, Partition_Index rank, Partition_Index procToWatch );
+void Test_FeEquationNumberRun_Regular( Dictionary* dictionary, void* context, IJK elSizes, 
+				       unsigned rank, unsigned nProcs, unsigned procToWatch );
 
 int main( int argc, char* argv[] ) {
-#if 0
 	MPI_Comm			CommWorld;
 	int				rank;
 	int				numProcessors;
@@ -156,7 +156,7 @@
 	MPI_Barrier( CommWorld );
 
 	elSizes[I_AXIS] = 3; elSizes[J_AXIS] = 3; elSizes[K_AXIS] = 3;
-	Test_FeEquationNumberRun_Regular( dictionary, context, "corner", elSizes, rank, procToWatch );
+	Test_FeEquationNumberRun_Regular( dictionary, context, elSizes, rank, numProcessors, procToWatch );
 
 	Stg_Class_Delete( context );
 	Stg_Class_Delete( dictionary );
@@ -172,104 +172,88 @@
 }
 
 
-void Test_FeEquationNumberRun_Regular( Dictionary* dictionary, void* context, const char* nLayoutType, IJK elSizes,
-	Partition_Index rank, Partition_Index procToWatch )
+FeEquationNumber* buildEqNum( unsigned nProcs, unsigned* sizes, Dictionary* dict, Context* context, 
+			      VariableCondition** bcs )
 {
-	Topology*			nTopology;
-	ElementLayout*			eLayout;
-	NodeLayout*			nLayout;
-	MeshDecomp*			decomp;
-	MeshLayout*			meshLayout;
-	ElementType_Register*		elementType_Register;
-	ExtensionManager_Register*		extensionMgr_Register;
-	FiniteElement_Mesh*		feMesh;
+	CartesianGenerator*		gen;
+	FeMesh*				feMesh;
+	DofLayout*			dofs;
+	FeEquationNumber*		eqNum;
+	Variable_Register*		varReg;
+	unsigned			maxDecomp[3] = {0, 1, 1};
+	double				minCrd[3];
+	double				maxCrd[3];
+	char*				varNames[] = {"x", "y", "z"};
+	LinkedDofInfo*			linkedDofInfo;
+
+	minCrd[0] = minCrd[1] = minCrd[2] = 0.0;
+	maxCrd[0] = maxCrd[1] = maxCrd[2] = 1.0;
+
+	gen = CartesianGenerator_New( "" );
+	CartesianGenerator_SetDimSize( gen, 3 );
+	CartesianGenerator_SetTopologyParams( gen, sizes, 0, NULL, maxDecomp );
+	CartesianGenerator_SetGeometryParams( gen, minCrd, maxCrd );
+	CartesianGenerator_SetShadowDepth( gen, 1 );
+
+	feMesh = FeMesh_New( "" );
+	Mesh_SetGenerator( feMesh, gen );
+	FeMesh_SetElementFamily( feMesh, "linear" );
+	Build( feMesh, NULL, False );
+
+	varReg = Variable_Register_New();
+
+	Variable_NewVector( "coords", Variable_DataType_Double, 3, 
+			    &feMesh->topo->domains[0]->nDomains, (void**)&feMesh->verts, 
+			    varReg, 
+			    varNames[0], 
+			    varNames[1], 
+			    varNames[2] );
+
+	dofs = DofLayout_New( "", varReg, 0, feMesh );
+	dofs->nBaseVariables = 3;
+	dofs->baseVariables = Memory_Alloc_Array_Unnamed( Variable*, 3 );
+	dofs->baseVariables[0] = Variable_Register_GetByName( varReg, varNames[0] );
+	dofs->baseVariables[1] = Variable_Register_GetByName( varReg, varNames[1] );
+	dofs->baseVariables[2] = Variable_Register_GetByName( varReg, varNames[2] );
+	Build( dofs, NULL, False );
+	Initialise( dofs, NULL, False );
+
+	*bcs = (VariableCondition*)WallVC_New( "WallVC", "wallBC", varReg, NULL, dict, feMesh );
+	Build( *bcs, NULL, False );
+	Initialise( *bcs, NULL, False );
+
+	linkedDofInfo = LinkedDofInfo_New( "linkedDofInfo", feMesh, dofs, dict );
+	Build( linkedDofInfo, context, False );
+
+	eqNum = FeEquationNumber_New( "feEquationNumber", feMesh, dofs, *bcs, linkedDofInfo );
+	Build( eqNum, NULL, False );
+	Initialise( eqNum, NULL, False );
+
+	return eqNum;
+}
+
+
+void Test_FeEquationNumberRun_Regular( Dictionary* dictionary, void* context, IJK elSizes, 
+				       unsigned rank, unsigned nProcs, unsigned procToWatch )
+{
 	VariableCondition*		vc;
 	FeEquationNumber*		feEquationNumber;
-	char*				varNames[] = {"x", "y", "z"};
-	Index				i, j;
-	Node_GlobalIndex		nodeCount;
-	Index				numVars = 3;
-	Variable_Register*		variableRegister;
-	DofLayout*			dofs;
 	Stream*				stream;
-	LinkedDofInfo*			linkedDofInfo = NULL;
-	int                             dim = Dictionary_GetUnsignedInt_WithDefault( dictionary, "dim", 3 );
-	#if 0
-	IndexSet*			bottomSet;
-	IndexSet*			leftSet;
-	IndexSet*			rightSet;
-	#endif
 
 	if( rank == procToWatch ) printf("Creating Geometry, Decomps and Layouts:\n");
 	Dictionary_Set( dictionary, "meshSizeI", Dictionary_Entry_Value_FromUnsignedInt( elSizes[0]+1 ) );
 	Dictionary_Set( dictionary, "meshSizeJ", Dictionary_Entry_Value_FromUnsignedInt( elSizes[1]+1 ) );
 	Dictionary_Set( dictionary, "meshSizeK", Dictionary_Entry_Value_FromUnsignedInt( elSizes[2]+1 ) );
-	/* create the layout, dof and mesh to use */
-	eLayout = (ElementLayout*)ParallelPipedHexaEL_New( "PPHexaEL", dim, dictionary );
-	if ("corner" == nLayoutType) {
-		nTopology = (Topology*)IJK6Topology_New( "IJK6Topology", dictionary );
-		nLayout = (NodeLayout*)CornerNL_New( "CornerNL", dictionary, eLayout, nTopology );
-	}
-	decomp = (MeshDecomp*)HexaMD_New( "HexaMD", dictionary, MPI_COMM_WORLD, eLayout, nLayout );
-	meshLayout = MeshLayout_New( "MeshLayout", eLayout, nLayout, decomp );
 
 	stream = Journal_Register (Info_Type, "myStream");
 
-	extensionMgr_Register = ExtensionManager_Register_New();
-	elementType_Register = ElementType_Register_New("elementTypeRegister");
-	ElementType_Register_Add( elementType_Register, (ElementType*)ConstantElementType_New("constant") );
-	ElementType_Register_Add( elementType_Register, (ElementType*)BilinearElementType_New("bilinear") );
-	ElementType_Register_Add( elementType_Register, (ElementType*)TrilinearElementType_New("trilinear") );
-	ElementType_Register_Add( elementType_Register, (ElementType*)LinearTriangleElementType_New("linear") );
 	if( rank == procToWatch ) printf("Creating F.E. Mesh:\n");
-	feMesh = FiniteElement_Mesh_New( "testMesh", meshLayout, sizeof(Node), sizeof(Element), extensionMgr_Register,
-		elementType_Register, dictionary );
-	
-	/* Create variable register */
-	variableRegister = Variable_Register_New();
-	
-	/* Create variables */
 	if( rank == procToWatch ) printf("Creating Vars:\n");
-	Variable_NewVector( 
-		"coords", 
-		Variable_DataType_Double, 
-		3, 
-		&feMesh->nodeLocalCount, 
-		(void**)&feMesh->nodeCoord, 
-		variableRegister, 
-		varNames[0], 
-		varNames[1], 
-		varNames[2] );
-
-
-	nodeCount = feMesh->layout->decomp->nodeLocalCount;
-	dofs = DofLayout_New( "dofLayout", variableRegister, nodeCount );
-
-	for (i = 0; i < nodeCount; i++)
-	{
-		for (j = 0; j < numVars; j++) {
-			DofLayout_AddDof_ByVarName(dofs, varNames[j], i);
-		}
-	}
-	Build(dofs, 0, False);
-	
 	if( rank == procToWatch ) printf("Creating VC:\n");
-	if ( IrregEL_Type == meshLayout->elementLayout->type) {
-		vc = (VariableCondition*) SetVC_New( "SetVC", "setBC", variableRegister, NULL, dictionary );
-	}
-	else {
-		vc = (VariableCondition*) WallVC_New( "WallVC", "wallBC", variableRegister, NULL, dictionary, feMesh );
-	}
-	
 	if( rank == procToWatch ) printf("Building:\n");
-	/* Build and initialise system */
 	if( rank == procToWatch ) printf("Building mesh:\n");
-	Build( feMesh, context, False );
 	if( rank == procToWatch ) printf("Building Variable Conditions:\n");
-	Build( vc, context, False );
-	
-	linkedDofInfo = LinkedDofInfo_New( "linkedDofInfo", feMesh, dofs, dictionary );
-	Build( linkedDofInfo, context, False );
+
 	/* Demonstration of setting up linked dof info through code */
 	#if 0
 	LinkedDofInfo_AddDofSet( linkedDofInfo );
@@ -284,45 +268,21 @@
 	#endif
 
 	if( rank == procToWatch ) printf("Creating EQ num:\n");
-	/* Create the finite element equation number utility */
-	feEquationNumber = FeEquationNumber_New( "feEquationNumber", feMesh, dofs, vc, linkedDofInfo );
-
-	
 	if( rank == procToWatch ) printf("Building FE Eq num:\n");
-	FeEquationNumber_Build( feEquationNumber );
-	
 	if( rank == procToWatch ) printf("Initialising:\n");
-	Initialise( feMesh, 0, False );
-	FeEquationNumber_Initialise( feEquationNumber );
-	
 	if( rank == procToWatch ) printf("Building LM:\n");
-	FeEquationNumber_BuildLocationMatrix( feEquationNumber );
+
+	feEquationNumber = buildEqNum( nProcs, elSizes, dictionary, context, &vc );
 	
 	if( rank == procToWatch ) {
-	
 		Journal_Printf( stream, "V.C. applied: " );
 		VariableCondition_PrintConcise( vc, stream );
 		FeEquationNumber_PrintDestinationArray( feEquationNumber, stream );
 		FeEquationNumber_PrintLocationMatrix( feEquationNumber, stream );
-		Print( linkedDofInfo, stream );
+		Print( feEquationNumber->linkedDofInfo, stream );
 	}
 	
-	
 	/* Destroy stuff */
 	if( rank == procToWatch ) printf("Cleaning Up:\n");
 	Stg_Class_Delete( feEquationNumber );
-	Stg_Class_Delete( vc );
-	Stg_Class_Delete( feMesh );
-	Stg_Class_Delete( extensionMgr_Register );
-	Stg_Class_Delete( elementType_Register );
-	Stg_Class_Delete( variableRegister );
-	Stg_Class_Delete( dofs );
-
-
-	Stg_Class_Delete( meshLayout );
-	Stg_Class_Delete( decomp );
-	Stg_Class_Delete( nLayout );
-	Stg_Class_Delete( eLayout );
-	Stg_Class_Delete( nTopology );
-#endif
 }



More information about the cig-commits mailing list