/[escript]/trunk/ripley/src/RipleyDomain.cpp
ViewVC logotype

Diff of /trunk/ripley/src/RipleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3755 by caltinay, Thu Jan 5 06:51:31 2012 UTC revision 3756 by caltinay, Fri Jan 6 02:35:19 2012 UTC
# Line 780  void RipleyDomain::updateTagsInUse(int f Line 780  void RipleyDomain::updateTagsInUse(int f
780      }      }
781  }  }
782    
783    //protected
784    Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
785            const IndexVector& index, const dim_t M, const dim_t N) const
786    {
787        // paso will manage the memory
788        index_t* indexC = MEMALLOC(index.size(), index_t);
789        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
790        copy(index.begin(), index.end(), indexC);
791        copy(ptr.begin(), ptr.end(), ptrC);
792        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
793    }
794    
795    //protected
796    Paso_Pattern* RipleyDomain::createMainPattern() const
797    {
798        IndexVector ptr(1,0);
799        IndexVector index;
800    
801        for (index_t i=0; i<getNumDOF(); i++) {
802            // add the DOF itself
803            index.push_back(i);
804            const dim_t num=insertNeighbourNodes(index, i);
805            ptr.push_back(ptr.back()+num+1);
806        }
807    
808        return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
809    }
810    
811    //protected
812    void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
813            const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
814    {
815        IndexVector ptr(1,0);
816        IndexVector index;
817        for (index_t i=0; i<getNumDOF(); i++) {
818            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
819            ptr.push_back(ptr.back()+colIndices[i].size());
820        }
821    
822        const dim_t M=ptr.size()-1;
823        *colPattern=createPasoPattern(ptr, index, M, N);
824    
825        IndexVector rowPtr(1,0);
826        IndexVector rowIndex;
827        for (dim_t id=0; id<N; id++) {
828            dim_t n=0;
829            for (dim_t i=0; i<M; i++) {
830                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
831                    if (index[j]==id) {
832                        rowIndex.push_back(i);
833                        n++;
834                        break;
835                    }
836                }
837            }
838            rowPtr.push_back(rowPtr.back()+n);
839        }
840    
841        // M and N are now reversed
842        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
843    }
844    
845    //protected
846    void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
847           const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
848           dim_t num_Sol, const vector<double>& array) const
849    {
850        const dim_t numMyCols = mat->pattern->mainPattern->numInput;
851        const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
852    
853        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
854        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
855        double* mainBlock_val = mat->mainBlock->val;
856        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
857        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
858        double* col_coupleBlock_val = mat->col_coupleBlock->val;
859        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
860        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
861        double* row_coupleBlock_val = mat->row_coupleBlock->val;
862    
863        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
864            // down columns of array
865            const dim_t j_Eq = nodes_Eq[k_Eq];
866            const dim_t i_row = j_Eq;
867    //printf("row:%d\n", i_row);
868            // only look at the matrix rows stored on this processor
869            if (i_row < numMyRows) {
870                for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
871                    const dim_t i_col = nodes_Sol[k_Sol];
872                    if (i_col < numMyCols) {
873                        for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
874                            if (mainBlock_index[k] == i_col) {
875                                mainBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
876                                break;
877                            }
878                        }
879                    } else {
880                        for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
881    //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;
882                            if (col_coupleBlock_index[k] == i_col - numMyCols) {
883                                col_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
884                                break;
885                            }
886                        }
887                    }
888                }
889            } else {
890                for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
891                    // across rows of array
892                    const dim_t i_col = nodes_Sol[k_Sol];
893                    if (i_col < numMyCols) {
894                        for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
895                             k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
896                        {
897    //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;
898                            if (row_coupleBlock_index[k] == i_col) {
899                                row_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
900                                break;
901                            }
902                        }
903                    }
904                }
905            }
906        }
907    }
908    
909  //  //
910  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
911  //  //
# Line 929  dim_t RipleyDomain::getNumDOF() const Line 1055  dim_t RipleyDomain::getNumDOF() const
1055      throw RipleyException("getNumDOF() not implemented");      throw RipleyException("getNumDOF() not implemented");
1056  }  }
1057    
1058    dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1059    {
1060        throw RipleyException("insertNeighbourNodes() not implemented");
1061    }
1062    
1063  void RipleyDomain::assembleCoordinates(escript::Data& arg) const  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1064  {  {
1065      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");

Legend:
Removed from v.3755  
changed lines
  Added in v.3756

  ViewVC Help
Powered by ViewVC 1.1.26