Index: ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp (revision 25485) +++ ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp (revision 25486) @@ -35,7 +35,7 @@ numberofdofspernode=femmodel->nodes->MaxNumDofs(GsetEnum); - /*if our matrices are coming from issm, we don't do dynamic allocation like Petsc + /*if our matrices are coming from issm, we don't do dynamic allocation like Petsc * does, and this routine is essentially useless. Force standard alloc in this case: */ toolkittype=ToolkitOptions::GetToolkitType(); @@ -117,6 +117,8 @@ int numnodesperelement = femmodel->elements->MaxNumNodes(); int numnodesperload = femmodel->loads->MaxNumNodes(); + int elementssize = femmodel->elements->Size(); + int loadssize = femmodel->loads->Size(); /*First, we are building chaining vectors so that we know what nodes are * connected to what elements. These vectors are such that: * for(int i=head[id];i!=-1;i=next[i]) @@ -123,11 +125,11 @@ * will loop over all the elements that are connected to the node number * id*/ head_e = xNew(localnumnodes); for(i=0;i(femmodel->elements->Size()*numnodesperelement); - count2offset_e = xNew(femmodel->elements->Size()*numnodesperelement); + next_e = xNew(elementssize*numnodesperelement); + count2offset_e = xNew(elementssize*numnodesperelement); k=0; - for(i=0;ielements->Size();i++){ + for(i=0;i(femmodel->elements->GetObjectByOffset(i)); int elementnumnodes = element->GetNumberOfNodes(); lidlist = xNew(elementnumnodes); @@ -148,10 +150,10 @@ /*Chain for loads*/ head_l = xNew(localnumnodes); for(i=0;i(femmodel->loads->Size()*numnodesperload); - count2offset_l = xNew(femmodel->loads->Size()*numnodesperload); + next_l = xNew(loadssize*numnodesperload); + count2offset_l = xNew(loadssize*numnodesperload); k=0; - for(i=0;iloads->Size();i++){ + for(i=0;i(femmodel->loads->GetObjectByOffset(i)); int loadnumnodes = load->GetNumberOfNodes(); lidlist = xNew(loadnumnodes); @@ -176,9 +178,11 @@ int *d_connectivity = xNewZeroInit(localnumnodes); int *o_connectivity = xNewZeroInit(localnumnodes); int flagsindices_counter; + int analysis_type; Vector *connectivity_clone= new Vector(localmasters,numnodes); + femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); /*Resetting flags to false at each iteration takes a lot of time, so we keep track of the flags * to reset in flagsindices, initialized with -1*/ for(i = 0;i(femmodel->nodes->GetObjectByOffset(i)); int lid = node->Lid(); int pid = node->Pid(); - /*Reinitialize flags to false*/ j=0; while(j(femmodel->elements->GetObjectByOffset(offset)); - element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,&flagsindices_counter,set1enum,set2enum); + element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,&flagsindices_counter,set1enum,set2enum,analysis_type); if(node->IsClone()){ connectivity_clone->SetValue(pid,d_nz+o_nz,ADD_VAL); } Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 25485) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 25486) @@ -736,6 +736,7 @@ Gauss* gauss=this->NewGauss(5); while(gauss->next()){ + this->JacobianDeterminant(&Jdet,xyz_list,gauss); /*Get strain rate assuming that epsilon has been allocated*/ @@ -3393,7 +3394,7 @@ } /*}}}*/ -void Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int* flagsindices_counter,int set1_enum,int set2_enum){/*{{{*/ +void Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int* flagsindices_counter,int set1_enum,int set2_enum, int analysis_type){/*{{{*/ /*Intermediaries*/ const int numnodes = this->GetNumberOfNodes(); @@ -3404,10 +3405,8 @@ /*Loop over all nodes*/ for(int i=0;inodes[i]->Lid(); if(!flags[nodelid]){ - /*flag current node so that no other element processes it*/ flags[nodelid]=true; @@ -3447,8 +3446,7 @@ /*Special case: 2d/3d coupling, the node of this element might be connected *to the basal element*/ - int analysis_type,approximation,numlayers; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); + int approximation,numlayers; if(analysis_type==StressbalanceAnalysisEnum){ this->GetInputValue(&approximation,ApproximationEnum); if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){ Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 25485) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 25486) @@ -163,7 +163,7 @@ void ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum); void ResultToMatrix(IssmDouble* values,int ncols,int output_enum); void ResultToVector(Vector* vector,int output_enum); - void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int* flagsindices_counter,int set1_enum,int set2_enum); + void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int* flagsindices_counter,int set1_enum,int set2_enum, int analysis_type); void SetBoolInput(Inputs* inputs,int enum_in,bool value); void SetIntInput(Inputs* inputs,int enum_in,int value);