/[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 3682 - (show annotations)
Mon Nov 21 01:56:45 2011 UTC (9 years, 8 months ago) by jfenwick
File size: 9009 byte(s)
Fixed some cases which allowed domain objects to be created without wrappers

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

  ViewVC Help
Powered by ViewVC 1.1.26