/[escript]/branches/refine/buckley/src/BuckleyDomain.cc
ViewVC logotype

Contents of /branches/refine/buckley/src/BuckleyDomain.cc

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3680 - (show annotations)
Fri Nov 18 04:48:53 2011 UTC (7 years, 8 months ago) by jfenwick
File size: 8551 byte(s)
It has memory errors during a getX.

1 #include "BuckleyDomain.h"
2 #include "BuckleyException.h"
3
4 #include <escript/FunctionSpace.h>
5
6 using namespace buckley;
7
8
9
10 namespace
11 {
12 // Need to put some thought into how many spaces we actually need
13 const int initialspace=1;
14 const int invalidspace=-1;
15
16 const int notIMPLEMENTED=-1;
17
18 const int ctsfn=initialspace;
19 }
20
21
22 BuckleyDomain::BuckleyDomain(double x, double y, double z)
23 :ot(x,y,z),modified(true),leaves(0)
24 {
25 }
26
27 BuckleyDomain::~BuckleyDomain()
28 {
29 if (leaves)
30 {
31 delete[] leaves;
32 }
33 }
34
35 bool BuckleyDomain::isValidFunctionSpaceType(int functionSpaceType) const
36 {
37 return functionSpaceType==initialspace;
38 }
39
40 escript::Data BuckleyDomain::getX() const
41 {
42 return escript::FunctionSpace(this->getPtr(), ctsfn).getX();
43 }
44
45
46 // typedef struct
47 // {
48 // escript::Vector* vec;
49 // size_t chunksize;
50 // unkid id;
51 //
52 // } idstruct;
53 //
54 //
55 // void copycoords(const OctCell& c, void* v)
56 // {
57 // idstruct* s=reinterpret_cast<idstruct*>(v);
58 // #ifdef OMP
59 // if (id / s->chunksize==omp_get_num_thread())
60 // {
61 // #endif
62 // s->vec[id++]
63 //
64 //
65 // #ifdef OMP
66 // }
67 // #endif
68 // }
69
70
71 void BuckleyDomain::setToX(escript::Data& arg) const
72 {
73 const BuckleyDomain& argDomain=dynamic_cast<const BuckleyDomain&>(*(arg.getFunctionSpace().getDomain()));
74 if (argDomain!=*this)
75 throw BuckleyException("Error - Illegal domain of data point locations");
76
77 if (arg.getFunctionSpace().getTypeCode()==invalidspace)
78 {
79 throw BuckleyException("Error - Invalid functionspace for coordinate loading");
80 }
81 if (modified) // is the cached data we have about this domain stale?
82 {
83 if (leaves!=0)
84 {
85 delete [] leaves;
86 }
87 leaves=ot.process();
88 }
89
90 if (arg.getFunctionSpace().getTypeCode()==getContinuousFunctionCode()) // values on nodes
91 {
92 escript::DataTypes::ValueType::ValueType vec=(&arg.getDataPointRW(0, 0));
93 int i; // one day OMP-3 will be default
94 // Why isn't this parallelised??
95 // well it isn't threadsafe the way it is written
96 // up to 8 cells could touch a point and try to write its coords
97 // This can be fixed but I haven't done it yet
98 // #pragma omp parallel for private(i) schedule(static)
99 for (i=0;i<ot.leafCount();++i)
100 {
101 const OctCell& oc=*leaves[i];
102 const double dx=oc.sides[0]/2;
103 const double dy=oc.sides[1]/2;
104 const double dz=oc.sides[2]/2;
105 const double cx=oc.centre[0];
106 const double cy=oc.centre[1];
107 const double cz=oc.centre[2];
108
109 // So this section is likely to be horribly slow
110 // Once the system is working unroll this
111 for (unsigned int j=0;j<8;++j)
112 {
113 unkid destpoint=oc.leafinfo->pmap[j];
114 if (destpoint<2)
115 {
116 continue; // it's a hanging node so don't produce
117 }
118
119 double x=cx;
120 double y=cy;
121 double z=cz;
122 if (j>3)
123 {
124 z+=dz;
125 }
126 else
127 {
128 z-=dz;
129 }
130 switch (j)
131 {
132 case 0:
133 case 3:
134 case 4:
135 case 7: x-=dx;
136
137 default: x+=dx;
138
139 }
140
141 switch (j)
142 {
143 case 0:
144 case 1:
145 case 4:
146 case 5: y-=dy;
147
148 default: y+=dy;
149 }
150
151 vec[destpoint*3]=x;
152 vec[destpoint*3+1]=y;
153 vec[destpoint*3+2]=z;
154 }
155
156 }
157 }
158 else // so it's not on the nodes and we need to do some interpolation
159 {
160 throw BuckleyException("Please don't use other functionspaces on this yet.");
161
162 }
163
164 // Dudley_Mesh* mesh=m_dudleyMesh.get();
165 // // in case of values node coordinates we can do the job directly:
166 // if (arg.getFunctionSpace().getTypeCode()==Nodes) {
167 // escriptDataC _arg=arg.getDataC();
168 // Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
169 // } else {
170 // escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
171 // escriptDataC _tmp_data=tmp_data.getDataC();
172 // Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
173 // // this is then interpolated onto arg:
174 // interpolateOnDomain(arg,tmp_data);
175 // }
176 // checkDudleyError();
177 }
178
179
180 std::string BuckleyDomain::getDescription() const
181 {
182 return "Dummy text for domain";
183 }
184
185 int BuckleyDomain::getContinuousFunctionCode() const
186 {
187 return ctsfn;
188 }
189
190 int BuckleyDomain::getReducedContinuousFunctionCode() const
191 {
192 return initialspace;
193 }
194
195 int BuckleyDomain::getFunctionCode() const
196 {
197 return invalidspace;
198 }
199
200 int BuckleyDomain::getReducedFunctionCode() const
201 {
202 return invalidspace;
203 }
204
205 int BuckleyDomain::getFunctionOnBoundaryCode() const
206 {
207 return invalidspace;
208 }
209
210 int BuckleyDomain::getReducedFunctionOnBoundaryCode() const
211 {
212 return invalidspace;
213 }
214
215 int BuckleyDomain::getFunctionOnContactZeroCode() const
216 {
217 return invalidspace;
218 }
219
220 int BuckleyDomain::getReducedFunctionOnContactZeroCode() const
221 {
222 return invalidspace;
223 }
224
225 int BuckleyDomain::getFunctionOnContactOneCode() const
226 {
227 return invalidspace;
228 }
229
230 int BuckleyDomain::getReducedFunctionOnContactOneCode() const
231 {
232 return invalidspace;
233 }
234
235 int BuckleyDomain::getSolutionCode() const
236 {
237 return initialspace;
238 }
239
240 int BuckleyDomain::getReducedSolutionCode() const
241 {
242 return initialspace;
243 }
244
245 // hahaha - no
246 int BuckleyDomain::getDiracDeltaFunctionsCode() const
247 {
248 return invalidspace ;
249 }
250
251 int BuckleyDomain::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
252 {
253 return notIMPLEMENTED;
254 }
255
256 int BuckleyDomain::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
257 {
258 return notIMPLEMENTED;
259 }
260
261 void BuckleyDomain::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
262 {
263 throw BuckleyException("Not Implemented");
264 }
265
266 void BuckleyDomain::addPDEToSystem(
267 escript::AbstractSystemMatrix& mat, escript::Data& rhs,
268 const escript::Data& A, const escript::Data& B, const escript::Data& C,
269 const escript::Data& D, const escript::Data& X, const escript::Data& Y,
270 const escript::Data& d, const escript::Data& y,
271 const escript::Data& d_contact, const escript::Data& y_contact,
272 const escript::Data& d_dirac, const escript::Data& y_dirac) const
273 {
274 throw BuckleyException("Not Implemented");
275 }
276
277 void BuckleyDomain::addPDEToRHS(escript::Data& rhs,
278 const escript::Data& X, const escript::Data& Y,
279 const escript::Data& y, const escript::Data& y_contact, const escript::Data& y_dirac) const
280 {
281 throw BuckleyException("Not Implemented");
282 }
283
284 void BuckleyDomain::addPDEToTransportProblem(
285 escript::AbstractTransportProblem& tp, escript::Data& source,
286 const escript::Data& M,
287 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,
288 const escript::Data& X,const escript::Data& Y,
289 const escript::Data& d, const escript::Data& y,
290 const escript::Data& d_contact,const escript::Data& y_contact,
291 const escript::Data& d_dirac,const escript::Data& y_dirac) const
292 {
293 throw BuckleyException("Not Implemented");
294 }
295
296 ASM_ptr BuckleyDomain::newSystemMatrix(
297 const int row_blocksize,
298 const escript::FunctionSpace& row_functionspace,
299 const int column_blocksize,
300 const escript::FunctionSpace& column_functionspace,
301 const int type) const
302 {
303 throw BuckleyException("Not Implemented");
304 }
305
306 ATP_ptr BuckleyDomain::newTransportProblem(
307 const bool useBackwardEuler,
308 const int blocksize,
309 const escript::FunctionSpace& functionspace,
310 const int type) const
311 {
312 throw BuckleyException("Not Implemented");
313 }
314
315 int BuckleyDomain::getNumDataPointsGlobal() const
316 {
317 throw BuckleyException("Not Implemented");
318 }
319
320
321 BUCKLEY_DLL_API
322 std::pair<int,int> BuckleyDomain::getDataShape(int functionSpaceCode) const
323 {
324 throw BuckleyException("Not Implemented");
325
326 }
327
328 BUCKLEY_DLL_API
329 void BuckleyDomain::setNewX(const escript::Data& arg)
330 {
331 throw BuckleyException("This domain does not support changing coordinates");
332 }
333
334
335 BUCKLEY_DLL_API
336 void BuckleyDomain::Print_Mesh_Info(const bool full) const
337 {
338 throw BuckleyException("Not Implemented");
339 }

  ViewVC Help
Powered by ViewVC 1.1.26