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

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

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

revision 3795 by caltinay, Thu Feb 2 05:54:34 2012 UTC revision 3800 by caltinay, Fri Feb 3 01:18:06 2012 UTC
# Line 649  void Rectangle::assembleGradient(escript Line 649  void Rectangle::assembleGradient(escript
649                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                      o[INDEX3(i,1,2,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
650                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                      o[INDEX3(i,0,3,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
651                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                      o[INDEX3(i,1,3,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
652                  } /* end of component loop i */                  } // end of component loop i
653              } /* end of k0 loop */              } // end of k0 loop
654          } /* end of k1 loop */          } // end of k1 loop
655      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedElements) {
656          out.requireWrite();          out.requireWrite();
657  #pragma omp parallel for  #pragma omp parallel for
# Line 665  void Rectangle::assembleGradient(escript Line 665  void Rectangle::assembleGradient(escript
665                  for (index_t i=0; i < numComp; ++i) {                  for (index_t i=0; i < numComp; ++i) {
666                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                      o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
667                      o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                      o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
668                  } /* end of component loop i */                  } // end of component loop i
669              } /* end of k0 loop */              } // end of k0 loop
670          } /* end of k1 loop */          } // end of k1 loop
671    
672      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == FaceElements) {
673          out.requireWrite();          out.requireWrite();
674  #pragma omp parallel  #pragma omp parallel
# Line 685  void Rectangle::assembleGradient(escript Line 686  void Rectangle::assembleGradient(escript
686                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
687                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
688                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
689                      } /* end of component loop i */                      } // end of component loop i
690                  } /* end of k1 loop */                  } // end of k1 loop
691              } /* end of face 0 */              } // end of face 0
692              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
693  #pragma omp for nowait  #pragma omp for nowait
694                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 701  void Rectangle::assembleGradient(escript Line 702  void Rectangle::assembleGradient(escript
702                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
703                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx3 + f_01[i]*cx1 + f_10[i]*cx4 + f_11[i]*cx6;
704                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,1,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
705                      } /* end of component loop i */                      } // end of component loop i
706                  } /* end of k1 loop */                  } // end of k1 loop
707              } /* end of face 1 */              } // end of face 1
708              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
709  #pragma omp for nowait  #pragma omp for nowait
710                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 717  void Rectangle::assembleGradient(escript Line 718  void Rectangle::assembleGradient(escript
718                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
719                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
720                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
721                      } /* end of component loop i */                      } // end of component loop i
722                  } /* end of k0 loop */                  } // end of k0 loop
723              } /* end of face 2 */              } // end of face 2
724              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
725  #pragma omp for nowait  #pragma omp for nowait
726                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 733  void Rectangle::assembleGradient(escript Line 734  void Rectangle::assembleGradient(escript
734                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy1 + f_01[i]*cy6 + f_10[i]*cy3 + f_11[i]*cy4;
735                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,1,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
736                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;                          o[INDEX3(i,1,1,numComp,2)] = f_00[i]*cy3 + f_01[i]*cy4 + f_10[i]*cy1 + f_11[i]*cy6;
737                      } /* end of component loop i */                      } // end of component loop i
738                  } /* end of k0 loop */                  } // end of k0 loop
739              } /* end of face 3 */              } // end of face 3
740          } // end of parallel section          } // end of parallel section
741    
742      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {      } else if (out.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
743          out.requireWrite();          out.requireWrite();
744  #pragma omp parallel  #pragma omp parallel
# Line 752  void Rectangle::assembleGradient(escript Line 754  void Rectangle::assembleGradient(escript
754                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
755                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
756                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_00[i]*cy0 + f_01[i]*cy7;
757                      } /* end of component loop i */                      } // end of component loop i
758                  } /* end of k1 loop */                  } // end of k1 loop
759              } /* end of face 0 */              } // end of face 0
760              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
761  #pragma omp for nowait  #pragma omp for nowait
762                  for (index_t k1=0; k1 < m_NE1; ++k1) {                  for (index_t k1=0; k1 < m_NE1; ++k1) {
# Line 766  void Rectangle::assembleGradient(escript Line 768  void Rectangle::assembleGradient(escript
768                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
769                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);                          o[INDEX3(i,0,0,numComp,2)] = cx5*(f_10[i] + f_11[i]) + cx2*(f_00[i] + f_01[i]);
770                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;                          o[INDEX3(i,1,0,numComp,2)] = f_10[i]*cy0 + f_11[i]*cy7;
771                      } /* end of component loop i */                      } // end of component loop i
772                  } /* end of k1 loop */                  } // end of k1 loop
773              } /* end of face 1 */              } // end of face 1
774              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
775  #pragma omp for nowait  #pragma omp for nowait
776                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 780  void Rectangle::assembleGradient(escript Line 782  void Rectangle::assembleGradient(escript
782                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
783                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_00[i]*cx0 + f_10[i]*cx7;
784                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);                          o[INDEX3(i,1,0,numComp,2)] = cy2*(f_00[i] + f_10[i]) + cy5*(f_01[i] + f_11[i]);
785                      } /* end of component loop i */                      } // end of component loop i
786                  } /* end of k0 loop */                  } // end of k0 loop
787              } /* end of face 2 */              } // end of face 2
788              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
789  #pragma omp for nowait  #pragma omp for nowait
790                  for (index_t k0=0; k0 < m_NE0; ++k0) {                  for (index_t k0=0; k0 < m_NE0; ++k0) {
# Line 794  void Rectangle::assembleGradient(escript Line 796  void Rectangle::assembleGradient(escript
796                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
797                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;                          o[INDEX3(i,0,0,numComp,2)] = f_01[i]*cx0 + f_11[i]*cx7;
798                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);                          o[INDEX3(i,1,0,numComp,2)] = cy5*(f_01[i] + f_11[i]) + cy2*(f_00[i] + f_10[i]);
799                      } /* end of component loop i */                      } // end of component loop i
800                  } /* end of k0 loop */                  } // end of k0 loop
801              } /* end of face 3 */              } // end of face 3
802          } // end of parallel section          } // end of parallel section
803      }      }
804  }  }
# Line 809  void Rectangle::assembleIntegrate(vector Line 811  void Rectangle::assembleIntegrate(vector
811      const double h1 = m_l1/m_gNE1;      const double h1 = m_l1/m_gNE1;
812      const index_t left = (m_offset0==0 ? 0 : 1);      const index_t left = (m_offset0==0 ? 0 : 1);
813      const index_t bottom = (m_offset1==0 ? 0 : 1);      const index_t bottom = (m_offset1==0 ? 0 : 1);
814      if (arg.getFunctionSpace().getTypeCode() == Elements) {      const int fs=arg.getFunctionSpace().getTypeCode();
815          const double w = h0*h1/4.;      if (fs == Elements && arg.actsExpanded()) {
816  #pragma omp parallel  #pragma omp parallel
817          {          {
818              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
819                const double w = h0*h1/4.;
820  #pragma omp for nowait  #pragma omp for nowait
821              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {              for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
822                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {                  for (index_t k0 = left; k0 < left+m_ownNE0; ++k0) {
# Line 824  void Rectangle::assembleIntegrate(vector Line 827  void Rectangle::assembleIntegrate(vector
827                          const double f2 = f[INDEX2(i,2,numComp)];                          const double f2 = f[INDEX2(i,2,numComp)];
828                          const double f3 = f[INDEX2(i,3,numComp)];                          const double f3 = f[INDEX2(i,3,numComp)];
829                          int_local[i]+=(f0+f1+f2+f3)*w;                          int_local[i]+=(f0+f1+f2+f3)*w;
830                      }  /* end of component loop i */                      }  // end of component loop i
831                  } /* end of k0 loop */                  } // end of k0 loop
832              } /* end of k1 loop */              } // end of k1 loop
   
833  #pragma omp critical  #pragma omp critical
834              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
835                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
836          } // end of parallel section          } // end of parallel section
837      } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {  
838        } else if (fs==ReducedElements || (fs==Elements && !arg.actsExpanded())) {
839          const double w = h0*h1;          const double w = h0*h1;
840  #pragma omp parallel  #pragma omp parallel
841          {          {
# Line 843  void Rectangle::assembleIntegrate(vector Line 846  void Rectangle::assembleIntegrate(vector
846                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));                      const double* f = arg.getSampleDataRO(INDEX2(k0, k1, m_NE0));
847                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
848                          int_local[i]+=f[i]*w;                          int_local[i]+=f[i]*w;
849                      }  /* end of component loop i */                      }
850                  } /* end of k0 loop */                  }
851              } /* end of k1 loop */              }
   
852  #pragma omp critical  #pragma omp critical
853              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
854                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
855          } // end of parallel section          } // end of parallel section
856      } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {  
857          const double w0 = h0/2.;      } else if (fs == FaceElements && arg.actsExpanded()) {
         const double w1 = h1/2.;  
858  #pragma omp parallel  #pragma omp parallel
859          {          {
860              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
861                const double w0 = h0/2.;
862                const double w1 = h1/2.;
863              if (m_faceOffset[0] > -1) {              if (m_faceOffset[0] > -1) {
864  #pragma omp for nowait  #pragma omp for nowait
865                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {                  for (index_t k1 = bottom; k1 < bottom+m_ownNE1; ++k1) {
# Line 865  void Rectangle::assembleIntegrate(vector Line 868  void Rectangle::assembleIntegrate(vector
868                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
869                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
870                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
871                      }  /* end of component loop i */                      }  // end of component loop i
872                  } /* end of k1 loop */                  } // end of k1 loop
873              }              }
874    
875              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
# Line 877  void Rectangle::assembleIntegrate(vector Line 880  void Rectangle::assembleIntegrate(vector
880                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
881                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
882                          int_local[i]+=(f0+f1)*w1;                          int_local[i]+=(f0+f1)*w1;
883                      }  /* end of component loop i */                      }  // end of component loop i
884                  } /* end of k1 loop */                  } // end of k1 loop
885              }              }
886    
887              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
# Line 889  void Rectangle::assembleIntegrate(vector Line 892  void Rectangle::assembleIntegrate(vector
892                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
893                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
894                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
895                      }  /* end of component loop i */                      }  // end of component loop i
896                  } /* end of k0 loop */                  } // end of k0 loop
897              }              }
898    
899              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
# Line 901  void Rectangle::assembleIntegrate(vector Line 904  void Rectangle::assembleIntegrate(vector
904                          const double f0 = f[INDEX2(i,0,numComp)];                          const double f0 = f[INDEX2(i,0,numComp)];
905                          const double f1 = f[INDEX2(i,1,numComp)];                          const double f1 = f[INDEX2(i,1,numComp)];
906                          int_local[i]+=(f0+f1)*w0;                          int_local[i]+=(f0+f1)*w0;
907                      }  /* end of component loop i */                      }  // end of component loop i
908                  } /* end of k0 loop */                  } // end of k0 loop
909              }              }
   
910  #pragma omp critical  #pragma omp critical
911              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
912                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
913          } // end of parallel section          } // end of parallel section
914      } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {  
915        } else if (fs==ReducedFaceElements || (fs==FaceElements && !arg.actsExpanded())) {
916  #pragma omp parallel  #pragma omp parallel
917          {          {
918              vector<double> int_local(numComp, 0);              vector<double> int_local(numComp, 0);
# Line 919  void Rectangle::assembleIntegrate(vector Line 922  void Rectangle::assembleIntegrate(vector
922                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[0]+k1);
923                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
924                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
925                      }  /* end of component loop i */                      }
926                  } /* end of k1 loop */                  }
927              }              }
928    
929              if (m_faceOffset[1] > -1) {              if (m_faceOffset[1] > -1) {
# Line 929  void Rectangle::assembleIntegrate(vector Line 932  void Rectangle::assembleIntegrate(vector
932                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);                      const double* f = arg.getSampleDataRO(m_faceOffset[1]+k1);
933                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
934                          int_local[i]+=f[i]*h1;                          int_local[i]+=f[i]*h1;
935                      }  /* end of component loop i */                      }
936                  } /* end of k1 loop */                  }
937              }              }
938    
939              if (m_faceOffset[2] > -1) {              if (m_faceOffset[2] > -1) {
# Line 939  void Rectangle::assembleIntegrate(vector Line 942  void Rectangle::assembleIntegrate(vector
942                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[2]+k0);
943                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
944                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
945                      }  /* end of component loop i */                      }
946                  } /* end of k0 loop */                  }
947              }              }
948    
949              if (m_faceOffset[3] > -1) {              if (m_faceOffset[3] > -1) {
# Line 949  void Rectangle::assembleIntegrate(vector Line 952  void Rectangle::assembleIntegrate(vector
952                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);                      const double* f = arg.getSampleDataRO(m_faceOffset[3]+k0);
953                      for (index_t i=0; i < numComp; ++i) {                      for (index_t i=0; i < numComp; ++i) {
954                          int_local[i]+=f[i]*h0;                          int_local[i]+=f[i]*h0;
955                      }  /* end of component loop i */                      }
956                  } /* end of k0 loop */                  }
957              }              }
958    
959  #pragma omp critical  #pragma omp critical
960              for (index_t i=0; i<numComp; i++)              for (index_t i=0; i<numComp; i++)
961                  integrals[i]+=int_local[i];                  integrals[i]+=int_local[i];
962          } // end of parallel section          } // end of parallel section
963      }      } // function space selector
964  }  }
965    
966  //protected  //protected

Legend:
Removed from v.3795  
changed lines
  Added in v.3800

  ViewVC Help
Powered by ViewVC 1.1.26