/[escript]/branches/trilinos_from_5897/dudley/src/Assemble_getAssembleParameters.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/Assemble_getAssembleParameters.cpp

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

revision 6008 by caltinay, Wed Feb 17 23:53:30 2016 UTC revision 6009 by caltinay, Wed Mar 2 04:13:26 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
17  /************************************************************************************/  /****************************************************************************
18    
19  /*    assemblage routines: prepares the assemble parameter set */    Assemblage routines: prepares the assemble parameter set
20    
21  /************************************************************************************/  *****************************************************************************/
   
 #define ESNEEDPYTHON  
 #include "esysUtils/first.h"  
22    
23  #include "Assemble.h"  #include "Assemble.h"
24  #include "ShapeTable.h"  #include "ShapeTable.h"
25    
26  #include <paso/SystemMatrix.h>  #include <paso/SystemMatrix.h>
27    
28  /************************************************************************************/  namespace dudley {
29    
30  void Dudley_Assemble_getAssembleParameters(Dudley_NodeFile* nodes,  void Assemble_getAssembleParameters(const Dudley_NodeFile* nodes,
31          Dudley_ElementFile* elements, escript::ASM_ptr mat,          const Dudley_ElementFile* elements, escript::ASM_ptr mat,
32          const escript::Data* F, bool reducedIntegrationOrder,          const escript::Data& F, bool reducedIntegrationOrder,
33          Dudley_Assemble_Parameters* parm)          Assemble_Parameters* parm)
34  {  {
35      paso::SystemMatrix* S = dynamic_cast<paso::SystemMatrix*>(mat.get());      paso::SystemMatrix* S = dynamic_cast<paso::SystemMatrix*>(mat.get());
36      if (mat && !S) {      if (mat && !S) {
37          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getAssembleParameters: Unknown matrix type.");          throw DudleyException("Assemble_getAssembleParameters: Unknown matrix type.");
         return;  
38      }      }
     Dudley_resetError();  
39      parm->shapeFns = NULL;      parm->shapeFns = NULL;
40      if (!isEmpty(F) && !isExpanded(F))      if (!F.isEmpty() && !F.actsExpanded()) {
41      {          throw DudleyException("Assemble_getAssembleParameters: Right hand side is not expanded.");
         Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getAssembleParameters: Right hand side is not expanded.");  
         return;  
42      }      }
43    
44      if (!getQuadShape(elements->numDim, reducedIntegrationOrder, &parm->shapeFns))      if (!getQuadShape(elements->numDim, reducedIntegrationOrder, &parm->shapeFns)) {
45      {          throw DudleyException("Assemble_getAssembleParameters: Can not locate shape functions.");
         Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getAssembleParameters: Can not locate shape functions.");  
46      }      }
47      /*  check the dimensions of S and F */      /*  check the dimensions of S and F */
48      if (S != NULL && !isEmpty(F))      if (S != NULL && !F.isEmpty()) {
49      {          if (!F.numSamplesEqual(1, (S->row_distribution->getMyNumComponents()*S->row_block_size) /
50          if (!numSamplesEqual               S->logical_row_block_size)) {
51              (F, 1,              throw DudleyException("Assemble_getAssembleParameters: number of rows of matrix and length of right hand side don't match.");
              (S->row_distribution->getMyNumComponents() * S->row_block_size) /  
              S->logical_row_block_size))  
         {  
             Dudley_setError(TYPE_ERROR,  
                             "Dudley_Assemble_getAssembleParameters: number of rows of matrix and length of right hand side don't match.");  
             return;  
52          }          }
53      }      }
54      /* get the number of equations and components */      /* get the number of equations and components */
55      if (S == NULL)      if (S == NULL) {
56      {          if (F.isEmpty()) {
         if (isEmpty(F))  
         {  
57              parm->numEqu = 1;              parm->numEqu = 1;
58              parm->numComp = 1;              parm->numComp = 1;
59          }          } else {
60          else              parm->numEqu = F.getDataPointSize();
         {  
             parm->numEqu = getDataPointSize(F);  
61              parm->numComp = parm->numEqu;              parm->numComp = parm->numEqu;
62          }          }
63      }      } else {
64      else          if (F.isEmpty()) {
     {  
         if (isEmpty(F))  
         {  
65              parm->numEqu = S->logical_row_block_size;              parm->numEqu = S->logical_row_block_size;
66              parm->numComp = S->logical_col_block_size;              parm->numComp = S->logical_col_block_size;
67          }          } else {
68          else              if (F.getDataPointSize() != S->logical_row_block_size) {
69          {                  throw DudleyException("Assemble_getAssembleParameters: matrix row block size and number of components of right hand side don't match.");
             if (getDataPointSize(F) != S->logical_row_block_size)  
             {  
                 Dudley_setError(TYPE_ERROR,  
                                 "Dudley_Assemble_getAssembleParameters: matrix row block size and number of components of right hand side don't match.");  
                 return;  
70              }              }
71              parm->numEqu = S->logical_row_block_size;              parm->numEqu = S->logical_row_block_size;
72              parm->numComp = S->logical_col_block_size;              parm->numComp = S->logical_col_block_size;
# Line 101  void Dudley_Assemble_getAssembleParamete Line 75  void Dudley_Assemble_getAssembleParamete
75      parm->col_DOF = nodes->degreesOfFreedomMapping->target;      parm->col_DOF = nodes->degreesOfFreedomMapping->target;
76      parm->row_DOF = nodes->degreesOfFreedomMapping->target;      parm->row_DOF = nodes->degreesOfFreedomMapping->target;
77      /* get the information for the labeling of the degrees of freedom from matrix */      /* get the information for the labeling of the degrees of freedom from matrix */
78      if (S != NULL)      if (S != NULL) {
79      {          // Make sure # rows in matrix == num DOF for one of:
80          /* Make sure # rows in matrix == num DOF for one of: full or reduced (use numLocalDOF for MPI) */          // full or reduced (use numLocalDOF for MPI)
81          if (S->row_distribution->getMyNumComponents() * S->row_block_size ==          if (S->row_distribution->getMyNumComponents() * S->row_block_size ==
82              parm->numEqu * nodes->degreesOfFreedomDistribution->getMyNumComponents())              parm->numEqu * nodes->degreesOfFreedomDistribution->getMyNumComponents()) {
         {  
83              parm->row_DOF_UpperBound = nodes->degreesOfFreedomDistribution->getMyNumComponents();              parm->row_DOF_UpperBound = nodes->degreesOfFreedomDistribution->getMyNumComponents();
84              parm->row_DOF = nodes->degreesOfFreedomMapping->target;              parm->row_DOF = nodes->degreesOfFreedomMapping->target;
85              parm->row_jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);              parm->row_jac = Dudley_ElementFile_borrowJacobians(elements, nodes, reducedIntegrationOrder);
86          }          } else if (S->row_distribution->getMyNumComponents() * S->row_block_size ==
87          else if (S->row_distribution->getMyNumComponents() * S->row_block_size ==                   parm->numEqu * nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents()) {
                  parm->numEqu * nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents())  
         {  
88              parm->row_DOF_UpperBound = nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents();              parm->row_DOF_UpperBound = nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents();
89              parm->row_DOF = nodes->reducedDegreesOfFreedomMapping->target;              parm->row_DOF = nodes->reducedDegreesOfFreedomMapping->target;
90              parm->row_jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);              parm->row_jac = Dudley_ElementFile_borrowJacobians(elements, nodes, reducedIntegrationOrder);
91          }          } else {
92          else              throw DudleyException("Assemble_getAssembleParameters: number of rows in matrix does not match the number of degrees of freedom in mesh");
         {  
             Dudley_setError(TYPE_ERROR,  
                             "Dudley_Assemble_getAssembleParameters: number of rows in matrix does not match the number of degrees of freedom in mesh");  
93          }          }
94          /* Make sure # cols in matrix == num DOF for one of: full or reduced (use numLocalDOF for MPI) */          // Make sure # cols in matrix == num DOF for one of:
95            // full or reduced (use numLocalDOF for MPI)
96          if (S->col_distribution->getMyNumComponents() * S->col_block_size ==          if (S->col_distribution->getMyNumComponents() * S->col_block_size ==
97              parm->numComp * nodes->degreesOfFreedomDistribution->getMyNumComponents())              parm->numComp * nodes->degreesOfFreedomDistribution->getMyNumComponents()) {
         {  
98              parm->col_DOF_UpperBound = nodes->degreesOfFreedomDistribution->getMyNumComponents();              parm->col_DOF_UpperBound = nodes->degreesOfFreedomDistribution->getMyNumComponents();
99              parm->col_DOF = nodes->degreesOfFreedomMapping->target;              parm->col_DOF = nodes->degreesOfFreedomMapping->target;
100              parm->row_jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);              parm->row_jac = Dudley_ElementFile_borrowJacobians(elements, nodes, reducedIntegrationOrder);
101          }          } else if (S->col_distribution->getMyNumComponents() * S->col_block_size ==
102          else if (S->col_distribution->getMyNumComponents() * S->col_block_size ==                   parm->numComp * nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents()) {
                  parm->numComp * nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents())  
         {  
103              parm->col_DOF_UpperBound = nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents();              parm->col_DOF_UpperBound = nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents();
104              parm->col_DOF = nodes->reducedDegreesOfFreedomMapping->target;              parm->col_DOF = nodes->reducedDegreesOfFreedomMapping->target;
105              parm->row_jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);              parm->row_jac = Dudley_ElementFile_borrowJacobians(elements, nodes, reducedIntegrationOrder);
106          }          } else {
107          else              throw DudleyException("Assemble_getAssembleParameters: number of columns in matrix does not match the number of degrees of freedom in mesh");
         {  
             Dudley_setError(TYPE_ERROR,  
                             "Dudley_Assemble_getAssembleParameters: number of columns in matrix does not match the number of degrees of freedom in mesh");  
108          }          }
109      }      }
110      if (!Dudley_noError())  
111          return;      // get the information from right hand side
112      /* get the information from right hand side */      if (!F.isEmpty()) {
113      if (!isEmpty(F))          if (F.numSamplesEqual(1, nodes->degreesOfFreedomDistribution->getMyNumComponents())) {
     {  
         if (numSamplesEqual(F, 1, nodes->degreesOfFreedomDistribution->getMyNumComponents()))  
         {  
114              parm->row_DOF_UpperBound = nodes->degreesOfFreedomDistribution->getMyNumComponents();              parm->row_DOF_UpperBound = nodes->degreesOfFreedomDistribution->getMyNumComponents();
115              parm->row_DOF = nodes->degreesOfFreedomMapping->target;              parm->row_DOF = nodes->degreesOfFreedomMapping->target;
116              parm->row_jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);              parm->row_jac = Dudley_ElementFile_borrowJacobians(elements, nodes, reducedIntegrationOrder);
117          }          } else if (F.numSamplesEqual(1, nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents())) {
         else if (numSamplesEqual  
                  (F, 1, nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents()))  
         {  
118              parm->row_DOF_UpperBound = nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents();              parm->row_DOF_UpperBound = nodes->reducedDegreesOfFreedomDistribution->getMyNumComponents();
119              parm->row_DOF = nodes->reducedDegreesOfFreedomMapping->target;              parm->row_DOF = nodes->reducedDegreesOfFreedomMapping->target;
120              parm->row_jac = Dudley_ElementFile_borrowJacobeans(elements, nodes, reducedIntegrationOrder);              parm->row_jac = Dudley_ElementFile_borrowJacobians(elements, nodes, reducedIntegrationOrder);
121            } else {
122                throw DudleyException("Assemble_getAssembleParameters: length of RHS vector does not match the number of degrees of freedom in mesh");
123          }          }
124          else          if (S == NULL) {
         {  
             Dudley_setError(TYPE_ERROR,  
                             "Dudley_Assemble_getAssembleParameters: length of RHS vector does not match the number of degrees of freedom in mesh");  
         }  
         if (S == NULL)  
         {  
125              parm->col_DOF_UpperBound = parm->row_DOF_UpperBound;              parm->col_DOF_UpperBound = parm->row_DOF_UpperBound;
126              parm->col_DOF = parm->row_DOF;              parm->col_DOF = parm->row_DOF;
127              parm->row_jac = parm->row_jac;              parm->row_jac = parm->row_jac;
128          }          }
129      }      }
130    
131      if (parm->row_jac->numDim != parm->row_jac->numDim)      if (parm->row_jac->numDim != parm->row_jac->numDim) {
132      {          throw DudleyException("Assemble_getAssembleParameters: spacial dimension for row and column shape function must match.");
         Dudley_setError(TYPE_ERROR,  
                         "Dudley_Assemble_getAssembleParameters: spacial dimension for row and column shape function must match.");  
133      }      }
134    
135      if (elements->numNodes < parm->row_jac->numShapes)      if (elements->numNodes < parm->row_jac->numShapes) {
136      {          throw DudleyException("Assemble_getAssembleParameters: too many nodes are expected by row.");
137          Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getAssembleParameters: too many nodes are expected by row.");      }
138      }      if (parm->row_jac->numElements != elements->numElements) {
139      if (parm->row_jac->numElements != elements->numElements)          throw DudleyException("Assemble_getAssembleParameters: number of elements for row is wrong.");
     {  
         Dudley_setError(TYPE_ERROR, "Dudley_Assemble_getAssembleParameters: number of elements for row is wrong.");  
140      }      }
141    
142      parm->numQuad = parm->row_jac->numQuad;      parm->numQuad = parm->row_jac->numQuad;
# Line 197  void Dudley_Assemble_getAssembleParamete Line 146  void Dudley_Assemble_getAssembleParamete
146      parm->numShapes = parm->row_jac->numShapes;      parm->numShapes = parm->row_jac->numShapes;
147  }  }
148    
149    } // namespace dudley
150    

Legend:
Removed from v.6008  
changed lines
  Added in v.6009

  ViewVC Help
Powered by ViewVC 1.1.26