Changeset 13881

Timestamp:

11/06/12 11:37:35 (12 years ago)

Author:

Mathieu Morlighem

Message:

NEW: finishing up stiffness matrix allocation algorithm

Location:

issm/trunk-jpl/src/c/classes

Files:

• FemModel.cpp (modified) (5 diffs)
• FemModel.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -340 +340 @@
 
         /*Intermediary*/
-        int      fsize,ssize;
-        int      connectivity, numberofdofspernode;
-        int      analysis_type, configuration_type;
+        int  fsize,ssize,flocalsize,slocalsize;
+        int  connectivity, numberofdofspernode;
+        int  analysis_type,configuration_type;
+        int  m,n,M,N;
+        int *d_nnz = NULL;
+        int *o_nnz = NULL;
 
         /*output*/
@@ -358 +361 @@
 
         /*retrieve node info*/
-        fsize=this->nodes->NumberOfDofs(configuration_type,FsetEnum);
-        ssize=this->nodes->NumberOfDofs(configuration_type,SsetEnum);
+        fsize      = this->nodes->NumberOfDofs(configuration_type,FsetEnum);
+        ssize      = this->nodes->NumberOfDofs(configuration_type,SsetEnum);
+        flocalsize = this->nodes->NumberOfDofsLocal(analysis_type,FsetEnum);
+        slocalsize = this->nodes->NumberOfDofsLocal(analysis_type,SsetEnum);
+
         numberofdofspernode=this->nodes->MaxNumDofs(configuration_type,GsetEnum);
 
@@ -369 +375 @@
         }
         else{
-                /*IN PROGRESS*/
+                /*Kff*/
+                m=flocalsize; n=flocalsize; /*local  sizes*/
+                M=fsize;      N=fsize;      /*global sizes*/
+                this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,FsetEnum);
+                Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
+                xDelete<int>(d_nnz);
+                xDelete<int>(o_nnz);
+
+                /*Kfs*/
+                m=flocalsize; n=slocalsize; /*local  sizes*/
+                M=fsize;      N=ssize;      /*global sizes*/
+                this->MatrixNonzeros(&d_nnz,&o_nnz,FsetEnum,SsetEnum);
+                Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
+                xDelete<int>(d_nnz);
+                xDelete<int>(o_nnz);
+
+                /*Vectors*/
+                df =new Vector<IssmDouble>(flocalsize,fsize);
+                pf =new Vector<IssmDouble>(flocalsize,fsize);
         }
 
@@ -379 +403 @@
 }
 /*}}}*/
+void FemModel::MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum){/*{{{*/
+
+        /*Intermediary*/
+        int      i,j,k,index,count;
+        int      analysis_type,configuration_type;
+        int      fsize,dim;
+        int      d_nz,o_nz;
+        Element *element = NULL;
+        int     *head    = NULL;
+        int     *next    = NULL;
+
+        /*output*/
+        int *d_nnz = NULL;
+        int *o_nnz = NULL;
+
+        /*retrive parameters: */
+        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        this->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        this->parameters->FindParam(&dim,MeshDimensionEnum);
+
+        /*Get vector size and number of nodes*/
+        int numnodes            = nodes->NumberOfNodes(analysis_type);
+        int numberofdofspernode = nodes->MaxNumDofs(configuration_type,GsetEnum);
+        int m                   = nodes->NumberOfDofsLocal(analysis_type,set1enum);
+
+        /*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])
+         * will loop over all the elements that are connected to the node number
+         * id*/
+        if(dim==2){
+                head=xNew<int>(numnodes); for(i=0;i<numnodes;i++) head[i]=-1;
+                next=xNew<int>(elements->Size()*3); /*3 = number of nodes per element*/
+        }
+        else if(dim==3){
+                head=xNew<int>(numnodes); for(i=0;i<numnodes;i++) head[i]=-1;
+                next=xNew<int>(elements->Size()*6); /*6 = number of nodes per element*/
+        }
+        else{
+                _error_("dim "<<dim<<" not supported yet");
+        }
+        k=0;
+        for(i=0;i<elements->Size();i++){
+                element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
+                for(int j=0;j<3;j++){
+                        int index =dynamic_cast<Tria*>(element)->nodes[j]->sid;//starts at 0 for a given analysis
+                        _assert_(k>=0 && k<numnodes*elements->Size() && index>=0 && index<numnodes);
+                        next[k]=head[index];
+                        head[index]=k++;
+                }
+        }
+
+        /*OK now count number of dofs and flag each nodes for each node i*/
+        bool *flags                  = xNew<bool>(numnodes);
+        int  *d_connectivity         = xNewZeroInit<int>(numnodes);
+        int  *o_connectivity         = xNewZeroInit<int>(numnodes);
+        int  *connectivity_clone     = xNewZeroInit<int>(numnodes);
+        int  *all_connectivity_clone = xNewZeroInit<int>(numnodes);
+
+        /*Create connectivity vector*/
+        for(i=0;i<nodes->Size();i++){
+                Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
+                if(node->InAnalysis(analysis_type)){
+
+                        /*Reinitialize flags to 0*/
+                        for(j=0;j<numnodes;j++) flags[j]=false;
+
+                        /*Loop over elements that hold node number i*/
+                        _assert_(head[node->Sid()]!=-1);
+                        for(j=head[node->Sid()];j!=-1;j=next[j]){
+                                if(dim==2){
+                                        index = (int)(double(j)/3);
+                                }
+                                else if(dim==3){
+                                        index = (int)(double(j)/6);
+                                }
+                                else{
+                                        _error_("dim "<<dim<<" not supported yet");
+                                }
+                                element=dynamic_cast<Element*>(elements->GetObjectByOffset(index));
+                                element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,set1enum,set2enum);
+                                if(node->IsClone()){
+                                        connectivity_clone[node->Sid()]+=d_nz+o_nz;
+                                }
+                                else{
+                                        d_connectivity[node->Sid()]+=d_nz;
+                                        o_connectivity[node->Sid()]+=o_nz;
+                                }
+                        }
+                }
+        }
+        xDelete<bool>(flags);
+        xDelete<int>(head);
+        xDelete<int>(next);
+
+        /*sum over all cpus*/
+        MPI_Allreduce((void*)connectivity_clone,(void*)all_connectivity_clone,numnodes,MPI_INT,MPI_SUM,IssmComm::GetComm());
+        xDelete<int>(connectivity_clone);
+
+        if(set1enum==FsetEnum){
+                count=0;
+                d_nnz=xNew<int>(m);
+                o_nnz=xNew<int>(m);
+                for(i=0;i<nodes->Size();i++){
+                        Node* node=dynamic_cast<Node*>(nodes->GetObjectByOffset(i));
+                        if(node->InAnalysis(analysis_type) && !node->IsClone()){
+                                for(j=0;j<node->indexing.fsize;j++){
+                                        _assert_(count<m);
+                                        d_nnz[count]=numberofdofspernode*(d_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
+                                        o_nnz[count]=numberofdofspernode*(o_connectivity[node->Sid()] + all_connectivity_clone[node->Sid()]);
+                                        if(d_nnz[count]>m)   d_nnz[count]=m;
+                                        if(o_nnz[count]>fsize-m) o_nnz[count]=fsize-m;
+                                        count++;
+                                }
+                        }
+                }
+                _assert_(m==count);
+        }
+        else{
+                _error_("STOP not implemented");
+        }
+        xDelete<int>(d_connectivity);
+        xDelete<int>(o_connectivity);
+        xDelete<int>(all_connectivity_clone);
+
+        /*Allocate ouptput pointer*/
+        *pd_nnz=d_nnz;
+        *po_nnz=o_nnz;
+
+}/*}}}*/
 int  FemModel::UpdateVertexPositionsx(void){ /*{{{*/
 
@@ -628 +782 @@
         Kfs->Assemble();
         pf->Assemble();
+        //Kff->AllocationInfo();
+        //Kfs->AllocationInfo();
+        //_error_("STOP");
 
         /*Assign output pointers: */

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -55 +55 @@
                 void Echo();
                 void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);
+                void MatrixNonzeros(int** pd_nnz,int** po_nnz,int set1enum,int set2enum);
                 void Solve(void);
                 void OutputResults(void);