/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Diff of /branches/diaplayground/ripley/src/Brick.cpp

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

revision 3711 by caltinay, Tue Dec 6 00:24:43 2011 UTC revision 3713 by caltinay, Tue Dec 6 04:43:29 2011 UTC
# Line 975  void Brick::setToGradient(escript::Data& Line 975  void Brick::setToGradient(escript::Data&
975      }      }
976  }  }
977    
978    void Brick::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
979    {
980        escript::Data& in = *const_cast<escript::Data*>(&arg);
981        const dim_t numComp = in.getDataPointSize();
982        const double h0 = m_l0/m_gNE0;
983        const double h1 = m_l1/m_gNE1;
984        const double h2 = m_l2/m_gNE2;
985        if (arg.getFunctionSpace().getTypeCode() == Elements) {
986            const double w_0 = h0*h1*h2/8.;
987    #pragma omp parallel
988            {
989                vector<double> int_local(numComp, 0);
990    #pragma omp for
991                for (index_t k2 = 0; k2 < m_NE2; ++k2) {
992                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
993                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
994                            const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
995                            for (index_t i=0; i < numComp; ++i) {
996                                const register double f_0 = f[INDEX2(i,0,numComp)];
997                                const register double f_1 = f[INDEX2(i,1,numComp)];
998                                const register double f_2 = f[INDEX2(i,2,numComp)];
999                                const register double f_3 = f[INDEX2(i,3,numComp)];
1000                                const register double f_4 = f[INDEX2(i,4,numComp)];
1001                                const register double f_5 = f[INDEX2(i,5,numComp)];
1002                                const register double f_6 = f[INDEX2(i,6,numComp)];
1003                                const register double f_7 = f[INDEX2(i,7,numComp)];
1004                                int_local[i]+=(f_0+f_1+f_2+f_3+f_4+f_5+f_6+f_7)*w_0;
1005                            }  /* end of component loop i */
1006                        } /* end of k0 loop */
1007                    } /* end of k1 loop */
1008                } /* end of k2 loop */
1009    
1010    #pragma omp critical
1011                for (index_t i=0; i<numComp; i++)
1012                    integrals[i]+=int_local[i];
1013            }
1014        } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
1015            const double w_0 = h0*h1*h2;
1016    #pragma omp parallel
1017            {
1018                vector<double> int_local(numComp, 0);
1019    #pragma omp for
1020                for (index_t k2 = 0; k2 < m_NE2; ++k2) {
1021                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
1022                        for (index_t k0 = 0; k0 < m_NE0; ++k0) {
1023                            const double* f = in.getSampleDataRO(INDEX3(k0, k1, k2, m_NE0, m_NE1));
1024                            for (index_t i=0; i < numComp; ++i) {
1025                                int_local[i]+=f[i]*w_0;
1026                            }  /* end of component loop i */
1027                        } /* end of k0 loop */
1028                    } /* end of k1 loop */
1029                } /* end of k2 loop */
1030    
1031    #pragma omp critical
1032                for (index_t i=0; i<numComp; i++)
1033                    integrals[i]+=int_local[i];
1034            }
1035        } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
1036            const double w_0 = h1*h2/4.;
1037            const double w_1 = h0*h2/4.;
1038            const double w_2 = h0*h1/4.;
1039    #pragma omp parallel
1040            {
1041                vector<double> int_local(numComp, 0);
1042                if (m_faceOffset[0] > -1) {
1043    #pragma omp for
1044                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1045                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1046                            const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1047                            for (index_t i=0; i < numComp; ++i) {
1048                                const register double f_0 = f[INDEX2(i,0,numComp)];
1049                                const register double f_1 = f[INDEX2(i,1,numComp)];
1050                                const register double f_2 = f[INDEX2(i,2,numComp)];
1051                                const register double f_3 = f[INDEX2(i,3,numComp)];
1052                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
1053                            }  /* end of component loop i */
1054                        } /* end of k1 loop */
1055                    } /* end of k2 loop */
1056                }
1057    
1058                if (m_faceOffset[1] > -1) {
1059    #pragma omp for
1060                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1061                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1062                            const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1063                            for (index_t i=0; i < numComp; ++i) {
1064                                const register double f_0 = f[INDEX2(i,0,numComp)];
1065                                const register double f_1 = f[INDEX2(i,1,numComp)];
1066                                const register double f_2 = f[INDEX2(i,2,numComp)];
1067                                const register double f_3 = f[INDEX2(i,3,numComp)];
1068                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
1069                            }  /* end of component loop i */
1070                        } /* end of k1 loop */
1071                    } /* end of k2 loop */
1072                }
1073    
1074                if (m_faceOffset[2] > -1) {
1075    #pragma omp for
1076                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1077                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1078                            const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1079                            for (index_t i=0; i < numComp; ++i) {
1080                                const register double f_0 = f[INDEX2(i,0,numComp)];
1081                                const register double f_1 = f[INDEX2(i,1,numComp)];
1082                                const register double f_2 = f[INDEX2(i,2,numComp)];
1083                                const register double f_3 = f[INDEX2(i,3,numComp)];
1084                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
1085                            }  /* end of component loop i */
1086                        } /* end of k1 loop */
1087                    } /* end of k2 loop */
1088                }
1089    
1090                if (m_faceOffset[3] > -1) {
1091    #pragma omp for
1092                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1093                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1094                            const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1095                            for (index_t i=0; i < numComp; ++i) {
1096                                const register double f_0 = f[INDEX2(i,0,numComp)];
1097                                const register double f_1 = f[INDEX2(i,1,numComp)];
1098                                const register double f_2 = f[INDEX2(i,2,numComp)];
1099                                const register double f_3 = f[INDEX2(i,3,numComp)];
1100                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_1;
1101                            }  /* end of component loop i */
1102                        } /* end of k1 loop */
1103                    } /* end of k2 loop */
1104                }
1105    
1106                if (m_faceOffset[4] > -1) {
1107    #pragma omp for
1108                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1109                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1110                            const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1111                            for (index_t i=0; i < numComp; ++i) {
1112                                const register double f_0 = f[INDEX2(i,0,numComp)];
1113                                const register double f_1 = f[INDEX2(i,1,numComp)];
1114                                const register double f_2 = f[INDEX2(i,2,numComp)];
1115                                const register double f_3 = f[INDEX2(i,3,numComp)];
1116                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
1117                            }  /* end of component loop i */
1118                        } /* end of k1 loop */
1119                    } /* end of k2 loop */
1120                }
1121    
1122                if (m_faceOffset[5] > -1) {
1123    #pragma omp for
1124                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1125                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1126                            const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1127                            for (index_t i=0; i < numComp; ++i) {
1128                                const register double f_0 = f[INDEX2(i,0,numComp)];
1129                                const register double f_1 = f[INDEX2(i,1,numComp)];
1130                                const register double f_2 = f[INDEX2(i,2,numComp)];
1131                                const register double f_3 = f[INDEX2(i,3,numComp)];
1132                                int_local[i]+=(f_0+f_1+f_2+f_3)*w_2;
1133                            }  /* end of component loop i */
1134                        } /* end of k1 loop */
1135                    } /* end of k2 loop */
1136                }
1137    
1138    #pragma omp critical
1139                for (index_t i=0; i<numComp; i++)
1140                    integrals[i]+=int_local[i];
1141            }
1142    
1143        } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
1144            const double w_0 = h1*h2;
1145            const double w_1 = h0*h2;
1146            const double w_2 = h0*h1;
1147    #pragma omp parallel
1148            {
1149                vector<double> int_local(numComp, 0);
1150                if (m_faceOffset[0] > -1) {
1151    #pragma omp for
1152                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1153                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1154                            const double* f = in.getSampleDataRO(m_faceOffset[0]+INDEX2(k1,k2,m_NE1));
1155                            for (index_t i=0; i < numComp; ++i) {
1156                                int_local[i]+=f[i]*w_0;
1157                            }  /* end of component loop i */
1158                        } /* end of k1 loop */
1159                    } /* end of k2 loop */
1160                }
1161    
1162                if (m_faceOffset[1] > -1) {
1163    #pragma omp for
1164                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1165                        for (index_t k1=0; k1 < m_NE1; ++k1) {
1166                            const double* f = in.getSampleDataRO(m_faceOffset[1]+INDEX2(k1,k2,m_NE1));
1167                            for (index_t i=0; i < numComp; ++i) {
1168                                int_local[i]+=f[i]*w_0;
1169                            }  /* end of component loop i */
1170                        } /* end of k1 loop */
1171                    } /* end of k2 loop */
1172                }
1173    
1174                if (m_faceOffset[2] > -1) {
1175    #pragma omp for
1176                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1177                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1178                            const double* f = in.getSampleDataRO(m_faceOffset[2]+INDEX2(k0,k2,m_NE0));
1179                            for (index_t i=0; i < numComp; ++i) {
1180                                int_local[i]+=f[i]*w_1;
1181                            }  /* end of component loop i */
1182                        } /* end of k1 loop */
1183                    } /* end of k2 loop */
1184                }
1185    
1186                if (m_faceOffset[3] > -1) {
1187    #pragma omp for
1188                    for (index_t k2=0; k2 < m_NE2; ++k2) {
1189                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1190                            const double* f = in.getSampleDataRO(m_faceOffset[3]+INDEX2(k0,k2,m_NE0));
1191                            for (index_t i=0; i < numComp; ++i) {
1192                                int_local[i]+=f[i]*w_1;
1193                            }  /* end of component loop i */
1194                        } /* end of k1 loop */
1195                    } /* end of k2 loop */
1196                }
1197    
1198                if (m_faceOffset[4] > -1) {
1199    #pragma omp for
1200                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1201                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1202                            const double* f = in.getSampleDataRO(m_faceOffset[4]+INDEX2(k0,k1,m_NE0));
1203                            for (index_t i=0; i < numComp; ++i) {
1204                                int_local[i]+=f[i]*w_2;
1205                            }  /* end of component loop i */
1206                        } /* end of k1 loop */
1207                    } /* end of k2 loop */
1208                }
1209    
1210                if (m_faceOffset[5] > -1) {
1211    #pragma omp for
1212                    for (index_t k1=0; k1 < m_NE1; ++k1) {
1213                        for (index_t k0=0; k0 < m_NE0; ++k0) {
1214                            const double* f = in.getSampleDataRO(m_faceOffset[5]+INDEX2(k0,k1,m_NE0));
1215                            for (index_t i=0; i < numComp; ++i) {
1216                                int_local[i]+=f[i]*w_2;
1217                            }  /* end of component loop i */
1218                        } /* end of k1 loop */
1219                    } /* end of k2 loop */
1220                }
1221    
1222    #pragma omp critical
1223                for (index_t i=0; i<numComp; i++)
1224                    integrals[i]+=int_local[i];
1225            }
1226    
1227        } else {
1228            stringstream msg;
1229            msg << "setToIntegrals() not implemented for "
1230                << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
1231            throw RipleyException(msg.str());
1232        }
1233    }
1234    
1235  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
1236                                              bool reducedColOrder) const                                              bool reducedColOrder) const
1237  {  {

Legend:
Removed from v.3711  
changed lines
  Added in v.3713

  ViewVC Help
Powered by ViewVC 1.1.26