00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #ifndef INTERPOLATOR_H
00012 #define INTERPOLATOR_H
00013
00014
00015 #include "Field/BareField.h"
00016 #include "Field/LField.h"
00017 #include "Field/CompressedBrickIterator.h"
00018 #include "Index/NDIndex.h"
00019 #include "AppTypes/Vektor.h"
00020 #include "Utility/IpplInfo.h"
00021
00022 #ifdef IPPL_USE_STANDARD_HEADERS
00023 #include <iostream>
00024 using namespace std;
00025 #else
00026 #include <iostream.h>
00027 #endif
00028
00029
00030
00031
00032 template <class C>
00033 class CenteringTag {
00034 #ifdef IPPL_PURIFY
00035
00036 public:
00037 CenteringTag() {}
00038 CenteringTag(const CenteringTag<C> &) {}
00039 CenteringTag<C>& operator=(const CenteringTag<C> &) { return *this; }
00040 #endif
00041 };
00042
00043
00044 template <class PT, unsigned Dim, class M>
00045 inline
00046 NDIndex<Dim> FindNGP(const M& mesh, const Vektor<PT,Dim>& ppos,
00047 CenteringTag<Cell>) {
00048 return mesh.getCellContaining(ppos);
00049 }
00050
00051 template <class PT, unsigned Dim, class M>
00052 inline
00053 NDIndex<Dim> FindNGP(const M& mesh, const Vektor<PT,Dim>& ppos,
00054 CenteringTag<Vert>) {
00055 return mesh.getNearestVertex(ppos);
00056 }
00057
00058
00059 template <class PT, unsigned Dim, class M>
00060 inline
00061 void FindPos(Vektor<PT,Dim>& pos, const M& mesh, const NDIndex<Dim>& indices,
00062 CenteringTag<Cell>) {
00063 pos = mesh.getCellPosition(indices);
00064 return;
00065 }
00066
00067 template <class PT, unsigned Dim, class M>
00068 inline
00069 void FindPos(Vektor<PT,Dim>& pos, const M& mesh, const NDIndex<Dim>& indices,
00070 CenteringTag<Vert>) {
00071 pos = mesh.getVertexPosition(indices);
00072 return;
00073 }
00074
00075
00076 template <class PT, unsigned Dim, class M>
00077 inline
00078 void FindDelta(Vektor<PT,Dim>& delta, const M& mesh, const NDIndex<Dim>& gp,
00079 CenteringTag<Cell>) {
00080 NDIndex<Dim> vp;
00081 for (unsigned d=0; d<Dim; ++d) vp[d] = gp[d] + 1;
00082 delta = mesh.getDeltaCell(vp);
00083 return;
00084 }
00085
00086 template <class PT, unsigned Dim, class M>
00087 inline
00088 void FindDelta(Vektor<PT,Dim>& delta, const M& mesh, const NDIndex<Dim>& gp,
00089 CenteringTag<Vert>) {
00090 delta = mesh.getDeltaVertex(gp);
00091 return;
00092 }
00093
00094
00095
00096
00097
00098
00099 template <class T, unsigned Dim, class InterpolatorType>
00100 struct InterpolatorTraits {};
00101
00102
00103 template <class T, unsigned Dim>
00104 struct CacheData1 {
00105 NDIndex<Dim> Index_m;
00106 Vektor<T,Dim> Delta_m;
00107 };
00108
00109 template <class T, unsigned Dim>
00110 ostream &operator<<(ostream &o, const CacheData1<T,Dim> &c)
00111 {
00112 o << "(" << c.Index_m << "," << c.Delta_m << ")";
00113 return o;
00114 }
00115
00116
00117
00118 template <class T, unsigned Dim>
00119 struct CacheDataCIC {
00120 NDIndex<Dim> Index_m;
00121 int Offset_m[Dim];
00122 Vektor<T,Dim> Delta_m;
00123 };
00124
00125 template <class T, unsigned Dim>
00126 ostream &operator<<(ostream &o, const CacheDataCIC<T,Dim> &c)
00127 {
00128 Vektor<int,Dim> offset;
00129 for (int i=0; i < Dim; ++i)
00130 offset[i] = c.Offset_m[i];
00131 o << "(" << c.Index_m << "," << c.Delta_m << "," << offset << ")";
00132 return o;
00133 }
00134
00135
00136
00137
00138
00139 class Interpolator {
00140
00141 protected:
00142
00143
00144 template <class T, unsigned Dim>
00145 static CompressedBrickIterator<T,Dim>
00146 getFieldIter(const BareField<T,Dim>& f, const NDIndex<Dim>& pt) {
00147 typename BareField<T,Dim>::const_iterator_if lf_i, lf_end = f.end_if();
00148 for (lf_i = f.begin_if(); lf_i != lf_end; ++lf_i) {
00149 LField<T,Dim>& lf(*(*lf_i).second);
00150 if ( lf.getOwned().contains(pt) ) {
00151
00152 return lf.begin(pt);
00153 }
00154 }
00155
00156
00157 for (lf_i = f.begin_if(); lf_i != lf_end; ++lf_i) {
00158 LField<T,Dim>& lf(*(*lf_i).second);
00159 if ( lf.getAllocated().contains(pt) ) {
00160
00161 return lf.begin(pt);
00162 }
00163 }
00164
00165
00166 ERRORMSG("Interpolator::getFieldIter: attempt to access non-local index");
00167 ERRORMSG(pt << " on node " << Ippl::myNode() << endl);
00168 ERRORMSG("Dumping local owned and allocated domains:" << endl);
00169 int lfc = 0;
00170 for ( lf_i = f.begin_if(); lf_i != lf_end ; ++lf_i, ++lfc ) {
00171 LField<T,Dim>& lf(*(*lf_i).second);
00172 ERRORMSG(lfc << ": owned = " << lf.getOwned());
00173 ERRORMSG(", allocated = " << lf.getAllocated() << endl);
00174 }
00175 ERRORMSG("Error occurred for BareField with layout = " << f.getLayout());
00176 ERRORMSG(endl);
00177 ERRORMSG("Calling abort ..." << endl);
00178 Ippl::abort();
00179 return (*(*(f.begin_if())).second).begin();
00180 }
00181
00182 public:
00183
00184 Interpolator() {}
00185 ~Interpolator() {}
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233 };
00234
00235 #endif // INTERPOLATOR_H
00236
00237
00238
00239
00240
00241
00242