/[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 3711 by caltinay, Tue Dec 6 00:24:43 2011 UTC revision 3713 by caltinay, Tue Dec 6 04:43:29 2011 UTC
# Line 510  void Rectangle::setToGradient(escript::D Line 510  void Rectangle::setToGradient(escript::D
510      }      }
511  }  }
512    
513    void Rectangle::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
514    {
515        escript::Data& in = *const_cast<escript::Data*>(&arg);
516        const dim_t numComp = in.getDataPointSize();
517        const double h0 = m_l0/m_gNE0;
518        const double h1 = m_l1/m_gNE1;
519        if (arg.getFunctionSpace().getTypeCode() == Elements) {
520            const double w_0 = h0*h1/4.;
521    #pragma omp parallel
522            {
523                vector<double> int_local(numComp, 0);
524    #pragma omp for
525                for (index_t k1 = 0; k1 < m_NE1; ++k1) {
526                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
527                        const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
528                        for (index_t i=0; i < numComp; ++i) {
529                            const register double f_0 = f[INDEX2(i,0,numComp)];
530                            const register double f_1 = f[INDEX2(i,1,numComp)];
531                            const register double f_2 = f[INDEX2(i,2,numComp)];
532                            const register double f_3 = f[INDEX2(i,3,numComp)];
533                            int_local[i]+=(f_0+f_1+f_2+f_3)*w_0;
534                        }  /* end of component loop i */
535                    } /* end of k0 loop */
536                } /* end of k1 loop */
537    
538    #pragma omp critical
539                for (index_t i=0; i<numComp; i++)
540                    integrals[i]+=int_local[i];
541            }
542        } else if (arg.getFunctionSpace().getTypeCode() == ReducedElements) {
543            const double w_0 = h0*h1;
544    #pragma omp parallel
545            {
546                vector<double> int_local(numComp, 0);
547    #pragma omp for
548                for (index_t k1 = 0; k1 < m_NE1; ++k1) {
549                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
550                        const double* f = in.getSampleDataRO(INDEX2(k0, k1, m_NE0));
551                        for (index_t i=0; i < numComp; ++i) {
552                            int_local[i]+=f[i]*w_0;
553                        }  /* end of component loop i */
554                    } /* end of k0 loop */
555                } /* end of k1 loop */
556    
557    #pragma omp critical
558                for (index_t i=0; i<numComp; i++)
559                    integrals[i]+=int_local[i];
560            }
561        } else if (arg.getFunctionSpace().getTypeCode() == FaceElements) {
562            const double w_0 = h0/2.;
563            const double w_1 = h1/2.;
564    #pragma omp parallel
565            {
566                vector<double> int_local(numComp, 0);
567                if (m_faceOffset[0] > -1) {
568    #pragma omp for
569                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
570                        const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
571                        for (index_t i=0; i < numComp; ++i) {
572                            const register double f_0 = f[INDEX2(i,0,numComp)];
573                            const register double f_1 = f[INDEX2(i,1,numComp)];
574                            int_local[i]+=(f_0+f_1)*w_1;
575                        }  /* end of component loop i */
576                    } /* end of k1 loop */
577                }
578    
579                if (m_faceOffset[1] > -1) {
580    #pragma omp for
581                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
582                        const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
583                        for (index_t i=0; i < numComp; ++i) {
584                            const register double f_0 = f[INDEX2(i,0,numComp)];
585                            const register double f_1 = f[INDEX2(i,1,numComp)];
586                            int_local[i]+=(f_0+f_1)*w_1;
587                        }  /* end of component loop i */
588                    } /* end of k1 loop */
589                }
590    
591                if (m_faceOffset[2] > -1) {
592    #pragma omp for
593                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
594                        const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
595                        for (index_t i=0; i < numComp; ++i) {
596                            const register double f_0 = f[INDEX2(i,0,numComp)];
597                            const register double f_1 = f[INDEX2(i,1,numComp)];
598                            int_local[i]+=(f_0+f_1)*w_0;
599                        }  /* end of component loop i */
600                    } /* end of k0 loop */
601                }
602    
603                if (m_faceOffset[3] > -1) {
604    #pragma omp for
605                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
606                        const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
607                        for (index_t i=0; i < numComp; ++i) {
608                            const register double f_0 = f[INDEX2(i,0,numComp)];
609                            const register double f_1 = f[INDEX2(i,1,numComp)];
610                            int_local[i]+=(f_0+f_1)*w_0;
611                        }  /* end of component loop i */
612                    } /* end of k0 loop */
613                }
614    
615    #pragma omp critical
616                for (index_t i=0; i<numComp; i++)
617                    integrals[i]+=int_local[i];
618            }
619        } else if (arg.getFunctionSpace().getTypeCode() == ReducedFaceElements) {
620    #pragma omp parallel
621            {
622                vector<double> int_local(numComp, 0);
623                if (m_faceOffset[0] > -1) {
624    #pragma omp for
625                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
626                        const double* f = in.getSampleDataRO(m_faceOffset[0]+k1);
627                        for (index_t i=0; i < numComp; ++i) {
628                            int_local[i]+=f[i]*h1;
629                        }  /* end of component loop i */
630                    } /* end of k1 loop */
631                }
632    
633                if (m_faceOffset[1] > -1) {
634    #pragma omp for
635                    for (index_t k1 = 0; k1 < m_NE1; ++k1) {
636                        const double* f = in.getSampleDataRO(m_faceOffset[1]+k1);
637                        for (index_t i=0; i < numComp; ++i) {
638                            int_local[i]+=f[i]*h1;
639                        }  /* end of component loop i */
640                    } /* end of k1 loop */
641                }
642    
643                if (m_faceOffset[2] > -1) {
644    #pragma omp for
645                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
646                        const double* f = in.getSampleDataRO(m_faceOffset[2]+k0);
647                        for (index_t i=0; i < numComp; ++i) {
648                            int_local[i]+=f[i]*h0;
649                        }  /* end of component loop i */
650                    } /* end of k0 loop */
651                }
652    
653                if (m_faceOffset[3] > -1) {
654    #pragma omp for
655                    for (index_t k0 = 0; k0 < m_NE0; ++k0) {
656                        const double* f = in.getSampleDataRO(m_faceOffset[3]+k0);
657                        for (index_t i=0; i < numComp; ++i) {
658                            int_local[i]+=f[i]*h0;
659                        }  /* end of component loop i */
660                    } /* end of k0 loop */
661                }
662    
663    #pragma omp critical
664                for (index_t i=0; i<numComp; i++)
665                    integrals[i]+=int_local[i];
666            }
667        } else {
668            stringstream msg;
669            msg << "setToIntegrals() not implemented for "
670                << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
671            throw RipleyException(msg.str());
672        }
673    }
674    
675  void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,  void Rectangle::addPDEToSystem(escript::AbstractSystemMatrix& mat,
676          escript::Data& rhs, const escript::Data& A, const escript::Data& B,          escript::Data& rhs, const escript::Data& A, const escript::Data& B,
677          const escript::Data& C, const escript::Data& D,          const escript::Data& C, const escript::Data& D,

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

  ViewVC Help
Powered by ViewVC 1.1.26