PLaSK library

Loading...
Searching...
No Matches
therm2d.cpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#include "therm2d.hpp"
15
16 namespace plask { namespace thermal { namespace tstatic {
17
18 template<typename Geometry2DType>
19 ThermalFem2DSolver<Geometry2DType>::ThermalFem2DSolver(const std::string& name) :
20 FemSolverWithMaskedMesh<Geometry2DType, RectangularMesh<2>>(name),
21 loopno(0),
22 outTemperature(this, &ThermalFem2DSolver<Geometry2DType>::getTemperatures),
23 outHeatFlux(this, &ThermalFem2DSolver<Geometry2DType>::getHeatFluxes),
24 outThermalConductivity(this, &ThermalFem2DSolver<Geometry2DType>::getThermalConductivity),
25 maxerr(0.05),
26 inittemp(300.)
27{
28 temperatures.reset();
29 fluxes.reset();
30 inHeat = 0.;
31}
32
33
34 template<typename Geometry2DType>
35 ThermalFem2DSolver<Geometry2DType>::~ThermalFem2DSolver() {
36}
37
38
39 template<typename Geometry2DType>
40 void ThermalFem2DSolver<Geometry2DType>::loadConfiguration(XMLReader &source, Manager &manager)
41{
42 while (source.requireTagOrEnd())
43 {
44 std::string param = source.getNodeName();
45
46 if (param == "temperature")
47 this->readBoundaryConditions(manager, source, temperature_boundary);
48
49 else if (param == "heatflux")
50 this->readBoundaryConditions(manager, source, heatflux_boundary);
51
52 else if (param == "convection")
53 this->readBoundaryConditions(manager, source, convection_boundary);
54
55 else if (param == "radiation")
56 this->readBoundaryConditions(manager, source, radiation_boundary);
57
58 else if (param == "loop") {
59 inittemp = source.getAttribute<double>("inittemp", inittemp);
60 maxerr = source.getAttribute<double>("maxerr", maxerr);
61 source.requireTagEnd();
62 }
63
64 else if (!this->parseFemConfiguration(source, manager)) {
65 this->parseStandardConfiguration(source, manager);
66 }
67 }
68}
69
70
71 template<typename Geometry2DType>
72 void ThermalFem2DSolver<Geometry2DType>::onInitialize() {
73 if (!this->geometry) throw NoGeometryException(this->getId());
74 if (!this->mesh) throw NoMeshException(this->getId());
75
76 FemSolverWithMaskedMesh<Geometry2DType, RectangularMesh<2>>::onInitialize();
77
78 loopno = 0;
79 temperatures.reset(this->maskedMesh->size(), inittemp);
80
81 thickness.reset(this->maskedMesh->getElementsCount(), NAN);
82 // Set stiffness matrix and load vector
83 for (auto elem: this->maskedMesh->elements())
84 {
85 if (!isnan(thickness[elem.getIndex()])) continue;
86 auto material = this->geometry->getMaterial(elem.getMidpoint());
87 double top = elem.getUpper1(), bottom = elem.getLower1();
88 size_t row = elem.getIndex1();
89 size_t itop = row+1, ibottom = row;
90 size_t c = elem.getIndex0();
91 for (size_t r = row; r > 0; r--) {
92 auto e = this->mesh->element(c, r-1);
93 auto m = this->geometry->getMaterial(e.getMidpoint());
94 if (m == material) { //TODO ignore doping
95 bottom = e.getLower1();
96 ibottom = r-1;
97 } else break;
98 }
99 for (size_t r = elem.getIndex1()+1; r < this->mesh->axis[1]->size()-1; r++) {
100 auto e = this->mesh->element(c, r);
101 auto m = this->geometry->getMaterial(e.getMidpoint());
102 if (m == material) { //TODO ignore doping
103 top = e.getUpper1();
104 itop = r+1;
105 } else break;
106 }
107 double h = top - bottom;
108 for (size_t r = ibottom; r != itop; ++r) {
109 size_t idx = this->maskedMesh->element(c, r).getIndex();
110 if (idx != RectangularMaskedMesh2D::Element::UNKNOWN_ELEMENT_INDEX)
111 thickness[idx] = h;
112 }
113 }
114}
115
116
117 template<typename Geometry2DType> void ThermalFem2DSolver<Geometry2DType>::onInvalidate() {
118 temperatures.reset();
119 fluxes.reset();
120 thickness.reset();
121}
122
123 enum BoundarySide { LEFT, RIGHT, TOP, BOTTOM };
124
138 template <typename ConditionT>
139 static void setBoundaries(const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,ConditionT>& boundary_conditions,
140 size_t i1, size_t i2, size_t i3, size_t i4, double width, double height,
141 double& F1, double& F2, double& F3, double& F4,
142 double& K11, double& K22, double& K33, double& K44,
143 double& K12, double& K23, double& K34, double& K41,
144 const std::function<double(double,ConditionT,ConditionT,size_t,size_t,BoundarySide)>& F_function,
145 const std::function<double(double,ConditionT,ConditionT,size_t,size_t,BoundarySide)>& Kmm_function,
146 const std::function<double(double,ConditionT,ConditionT,size_t,size_t,BoundarySide)>& Kmn_function
147 )
148{
149 auto val1 = boundary_conditions.getValue(i1);
150 auto val2 = boundary_conditions.getValue(i2);
151 auto val3 = boundary_conditions.getValue(i3);
152 auto val4 = boundary_conditions.getValue(i4);
153 if (val1 && val2) { // bottom
154 F1 += F_function(width, *val1, *val2, i1, i2, BOTTOM); F2 += F_function(width, *val2, *val1, i2, i1, BOTTOM);
155 K11 += Kmm_function(width, *val1, *val2, i1, i2, BOTTOM); K22 += Kmm_function(width, *val2, *val1, i2, i1, BOTTOM);
156 K12 += Kmn_function(width, *val1, *val2, i1, i2, BOTTOM);
157 }
158 if (val2 && val3) { // right
159 F2 += F_function(height, *val2, *val3, i2, i3, RIGHT); F3 += F_function(height, *val3, *val2, i3, i2, RIGHT);
160 K22 += Kmm_function(height, *val2, *val3, i2, i3, RIGHT); K33 += Kmm_function(height, *val3, *val2, i3, i2, RIGHT);
161 K23 += Kmn_function(height, *val2, *val3, i2, i3, RIGHT);
162 }
163 if (val3 && val4) { // top
164 F3 += F_function(width, *val3, *val4, i3, i4, TOP); F4 += F_function(width, *val4, *val3, i4, i3, TOP);
165 K33 += Kmm_function(width, *val3, *val4, i3, i4, TOP); K44 += Kmm_function(width, *val4, *val3, i4, i3, TOP);
166 K34 += Kmn_function(width, *val3, *val4, i3, i4, TOP);
167 }
168 if (val4 && val1) { // left
169 F1 += F_function(height, *val1, *val4, i1, i4, LEFT); F4 += F_function(height, *val4, *val1, i4, i1, LEFT);
170 K11 += Kmm_function(height, *val1, *val4, i1, i4, LEFT); K44 += Kmm_function(height, *val4, *val1, i4, i1, LEFT);
171 K41 += Kmn_function(height, *val1, *val4, i1, i4, LEFT);
172 }
173}
174
175
176 template<>
177 void ThermalFem2DSolver<Geometry2DCartesian>::setMatrix(FemMatrix& A, DataVector<double>& B,
178 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,double>& btemperature,
179 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,double>& bheatflux,
180 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,Convection>& bconvection,
181 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,Radiation>& bradiation
182 )
183{
184 this->writelog(LOG_DETAIL, "Setting up matrix system ({})", A.describe());
185
186 auto iMesh = (this->maskedMesh)->getElementMesh();
187 auto heatdensities = inHeat(iMesh);
188
189 A.clear();
190 B.fill(0.);
191
192 // Set stiffness matrix and load vector
193 for (auto elem: this->maskedMesh->elements())
194 {
195 // nodes numbers for the current element
196 size_t loleftno = elem.getLoLoIndex();
197 size_t lorghtno = elem.getUpLoIndex();
198 size_t upleftno = elem.getLoUpIndex();
199 size_t uprghtno = elem.getUpUpIndex();
200
201 // element size
202 double elemwidth = elem.getUpper0() - elem.getLower0();
203 double elemheight = elem.getUpper1() - elem.getLower1();
204
205 // point and material in the middle of the element
206 Vec<2,double> midpoint = elem.getMidpoint();
207 auto material = this->geometry->getMaterial(midpoint);
208
209 // average temperature on the element
210 double temp = 0.25 * (temperatures[loleftno] + temperatures[lorghtno] + temperatures[upleftno] + temperatures[uprghtno]);
211
212 // thermal conductivity
213 double kx, ky;
214 std::tie(kx,ky) = std::tuple<double,double>(material->thermk(temp, thickness[elem.getIndex()]));
215
216 kx *= elemheight; kx /= elemwidth;
217 ky *= elemwidth; ky /= elemheight;
218
219 // load vector: heat densities
220 double f = 0.25e-12 * elemwidth * elemheight * heatdensities[elem.getIndex()]; // 1e-12 -> to transform µm² into m²
221
222 // set symmetric matrix components
223 double k44, k33, k22, k11, k43, k21, k42, k31, k32, k41;
224
225 k44 = k33 = k22 = k11 = (kx + ky) / 3.;
226 k43 = k21 = (-2. * kx + ky) / 6.;
227 k42 = k31 = - (kx + ky) / 6.;
228 k32 = k41 = (kx - 2. * ky) / 6.;
229
230 double f1 = f, f2 = f, f3 = f, f4 = f;
231
232 // boundary conditions: heat flux
233 setBoundaries<double>(bheatflux, loleftno, lorghtno, uprghtno, upleftno, elemwidth, elemheight,
234 f1, f2, f3, f4, k11, k22, k33, k44, k21, k32, k43, k41,
235 [](double len, double val, double, size_t, size_t, BoundarySide) { // F
236 return - 0.5e-6 * len * val;
237 },
238 [](double, double, double, size_t, size_t, BoundarySide){return 0.;}, // K diagonal
239 [](double, double, double, size_t, size_t, BoundarySide){return 0.;} // K off-diagonal
240 );
241
242 // boundary conditions: convection
243 setBoundaries<Convection>(bconvection, loleftno, lorghtno, uprghtno, upleftno, elemwidth, elemheight,
244 f1, f2, f3, f4, k11, k22, k33, k44, k21, k32, k43, k41,
245 [](double len, Convection val, Convection, size_t, size_t, BoundarySide) { // F
246 return 0.5e-6 * len * val.coeff * val.ambient;
247 },
248 [](double len, Convection val1, Convection val2, size_t, size_t, BoundarySide) { // K diagonal
249 return (val1.coeff + val2.coeff) * len / 6.;
250 },
251 [](double len, Convection val1, Convection val2, size_t, size_t, BoundarySide) { // K off-diagonal
252 return (val1.coeff + val2.coeff) * len / 12.;
253 }
254 );
255
256 // boundary conditions: radiation
257 setBoundaries<Radiation>(bradiation, loleftno, lorghtno, uprghtno, upleftno, elemwidth, elemheight,
258 f1, f2, f3, f4, k11, k22, k33, k44, k21, k32, k43, k41,
259 [this](double len, Radiation val, Radiation, size_t i, size_t, BoundarySide) -> double { // F
260 double a = val.ambient; a = a*a;
261 double T = this->temperatures[i]; T = T*T;
262 return - 0.5e-6 * len * val.emissivity * phys::SB * (T*T - a*a);},
263 [](double, Radiation, Radiation, size_t, size_t, BoundarySide){return 0.;}, // K diagonal
264 [](double, Radiation, Radiation, size_t, size_t, BoundarySide){return 0.;} // K off-diagonal
265 );
266
267 // set stiffness matrix
268 A(loleftno, loleftno) += k11;
269 A(lorghtno, lorghtno) += k22;
270 A(uprghtno, uprghtno) += k33;
271 A(upleftno, upleftno) += k44;
272
273 A(lorghtno, loleftno) += k21;
274 A(uprghtno, loleftno) += k31;
275 A(upleftno, loleftno) += k41;
276 A(uprghtno, lorghtno) += k32;
277 A(upleftno, lorghtno) += k42;
278 A(upleftno, uprghtno) += k43;
279
280 // set load vector
281 B[loleftno] += f1;
282 B[lorghtno] += f2;
283 B[uprghtno] += f3;
284 B[upleftno] += f4;
285 }
286
287 // boundary conditions of the first kind
288 A.applyBC(btemperature, B);
289
290#ifndef NDEBUG
291 double* aend = A.data + A.size;
292 for (double* pa = A.data; pa != aend; ++pa) {
293 if (isnan(*pa) || isinf(*pa))
294 throw ComputationError(this->getId(), "error in stiffness matrix at position {0}", pa-A.data);
295 }
296#endif
297
298}
299
300
301 template<>
302 void ThermalFem2DSolver<Geometry2DCylindrical>::setMatrix(FemMatrix& A, DataVector<double>& B,
303 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,double>& btemperature,
304 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,double>& bheatflux,
305 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,Convection>& bconvection,
306 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,Radiation>& bradiation
307 )
308{
309 this->writelog(LOG_DETAIL, "Setting up matrix system ({})", A.describe());
310
311 auto iMesh = (this->maskedMesh)->getElementMesh();
312 auto heatdensities = inHeat(iMesh);
313
314 A.clear();
315 B.fill(0.);
316
317 // Set stiffness matrix and load vector
318 for (auto elem: this->maskedMesh->elements())
319 {
320 // nodes numbers for the current element
321 size_t loleftno = elem.getLoLoIndex();
322 size_t lorghtno = elem.getUpLoIndex();
323 size_t upleftno = elem.getLoUpIndex();
324 size_t uprghtno = elem.getUpUpIndex();
325
326 // element size
327 double elemwidth = elem.getUpper0() - elem.getLower0();
328 double elemheight = elem.getUpper1() - elem.getLower1();
329
330 // point and material in the middle of the element
331 Vec<2,double> midpoint = elem.getMidpoint();
332 auto material = geometry->getMaterial(midpoint);
333 double r = midpoint.rad_r();
334
335 // average temperature on the element
336 double temp = 0.25 * (temperatures[loleftno] + temperatures[lorghtno] + temperatures[upleftno] + temperatures[uprghtno]);
337
338 // thermal conductivity
339 double kx, ky;
340 std::tie(kx,ky) = std::tuple<double,double>(material->thermk(temp, thickness[elem.getIndex()]));
341
342 kx = kx * elemheight / elemwidth;
343 ky = ky * elemwidth / elemheight;
344
345 // load vector: heat densities
346 double f = 0.25e-12 * r * elemwidth * elemheight * heatdensities[elem.getIndex()]; // 1e-12 -> to transform µm² into m²
347
348 // set symmetric matrix components
349 double k44, k33, k22, k11, k43, k21, k42, k31, k32, k41;
350
351 k44 = k33 = k22 = k11 = (kx + ky) / 3.;
352 k43 = k21 = (-2. * kx + ky) / 6.;
353 k42 = k31 = - (kx + ky) / 6.;
354 k32 = k41 = (kx - 2. * ky) / 6.;
355
356 double f1 = f, f2 = f, f3 = f, f4 = f;
357
358 // boundary conditions: heat flux
359 setBoundaries<double>(bheatflux, loleftno, lorghtno, uprghtno, upleftno, elemwidth, elemheight,
360 f1, f2, f3, f4, k11, k22, k33, k44, k21, k32, k43, k41,
361 [&](double len, double val, double, size_t i1, size_t i2, BoundarySide side) -> double { // F
362 if (side == LEFT) return - 0.5e-6 * len * val * elem.getLower0();
363 else if (side == RIGHT) return - 0.5e-6 * len * val * elem.getUpper0();
364 else return - 0.5e-6 * len * val * (r + (i1<i2? -len/6. : len/6.));
365 },
366 [](double, double, double, size_t, size_t, BoundarySide){return 0.;}, // K diagonal
367 [](double, double, double, size_t, size_t, BoundarySide){return 0.;} // K off-diagonal
368 );
369
370 // boundary conditions: convection
371 setBoundaries<Convection>(bconvection, loleftno, lorghtno, uprghtno, upleftno, elemwidth, elemheight,
372 f1, f2, f3, f4, k11, k22, k33, k44, k21, k32, k43, k41,
373 [&](double len, Convection val1, Convection val2, size_t i1, size_t i2, BoundarySide side) -> double { // F
374 double a = 0.125e-6 * len * (val1.coeff + val2.coeff) * (val1.ambient + val2.ambient);
375 if (side == LEFT) return a * elem.getLower0();
376 else if (side == RIGHT) return a * elem.getUpper0();
377 else return a * (r + (i1<i2? -len/6. : len/6.));
378
379 },
380 [&](double len, Convection val1, Convection val2, size_t i1, size_t i2, BoundarySide side) -> double { // K diagonal
381 double a = (val1.coeff + val2.coeff) * len / 6.;
382 if (side == LEFT) return a * elem.getLower0();
383 else if (side == RIGHT) return a * elem.getUpper0();
384 else return a * (r + (i1<i2? -len/6. : len/6.));
385 },
386 [&](double len, Convection val1, Convection val2, size_t, size_t, BoundarySide side) -> double { // K off-diagonal
387 double a = (val1.coeff + val2.coeff) * len / 12.;
388 if (side == LEFT) return a * elem.getLower0();
389 else if (side == RIGHT) return a * elem.getUpper0();
390 else return a * r;
391 }
392 );
393
394 // boundary conditions: radiation
395 setBoundaries<Radiation>(bradiation, loleftno, lorghtno, uprghtno, upleftno, elemwidth, elemheight,
396 f1, f2, f3, f4, k11, k22, k33, k44, k21, k32, k43, k41,
397 [&,this](double len, Radiation val, Radiation, size_t i1, size_t i2, BoundarySide side) -> double { // F
398 double amb = val.ambient; amb = amb*amb;
399 double T = this->temperatures[i1]; T = T*T;
400 double a = - 0.5e-6 * len * val.emissivity * phys::SB * (T*T - amb*amb);
401 if (side == LEFT) return a * elem.getLower0();
402 else if (side == RIGHT) return a * elem.getUpper0();
403 else return a * (r + (i1<i2? -len/6. : len/6.));
404 },
405 [](double, Radiation, Radiation, size_t, size_t, BoundarySide){return 0.;}, // K diagonal
406 [](double, Radiation, Radiation, size_t, size_t, BoundarySide){return 0.;} // K off-diagonal
407 );
408
409 // set stiffness matrix
410 A(loleftno, loleftno) += r * k11;
411 A(lorghtno, lorghtno) += r * k22;
412 A(uprghtno, uprghtno) += r * k33;
413 A(upleftno, upleftno) += r * k44;
414
415 A(lorghtno, loleftno) += r * k21;
416 A(uprghtno, loleftno) += r * k31;
417 A(upleftno, loleftno) += r * k41;
418 A(uprghtno, lorghtno) += r * k32;
419 A(upleftno, lorghtno) += r * k42;
420 A(upleftno, uprghtno) += r * k43;
421
422 // set load vector
423 B[loleftno] += f1;
424 B[lorghtno] += f2;
425 B[uprghtno] += f3;
426 B[upleftno] += f4;
427 }
428
429 // boundary conditions of the first kind
430 A.applyBC(btemperature, B);
431
432#ifndef NDEBUG
433 double* aend = A.data + A.size;
434 for (double* pa = A.data; pa != aend; ++pa) {
435 if (isnan(*pa) || isinf(*pa))
436 throw ComputationError(this->getId(), "error in stiffness matrix at position {0}", pa-A.data);
437 }
438#endif
439
440}
441
442
443 template<typename Geometry2DType>
444 double ThermalFem2DSolver<Geometry2DType>::compute(int loops){
445 this->initCalculation();
446
447 fluxes.reset();
448
449 // store boundary conditions for current mesh
450 auto btemperature = temperature_boundary(this->maskedMesh, this->geometry);
451 auto bheatflux = heatflux_boundary(this->maskedMesh, this->geometry);
452 auto bconvection = convection_boundary(this->maskedMesh, this->geometry);
453 auto bradiation = radiation_boundary(this->maskedMesh, this->geometry);
454
455 this->writelog(LOG_INFO, "Running thermal calculations");
456
457 int loop = 0;
458 size_t size = this->maskedMesh->size();
459
460 std::unique_ptr<FemMatrix> pA(this->getMatrix());
461 FemMatrix& A = *pA.get();
462
463 double err;
464 toterr = 0.;
465
466# ifndef NDEBUG
467 if (!temperatures.unique()) this->writelog(LOG_DEBUG, "Temperature data held by something else...");
468# endif
469 temperatures = temperatures.claim();
470
471 DataVector<double> B(size), temp0(size);
472
473 do {
474 std::copy(temperatures.begin(), temperatures.end(), temp0.begin());
475
476 setMatrix(A, B, btemperature, bheatflux, bconvection, bradiation);
477 A.solve(B, temperatures);
478
479 // Update error
480 err = 0.;
481 maxT = 0.;
482 for (auto temp = temperatures.begin(), t = temp0.begin(); t != temp0.end(); ++temp, ++t)
483 {
484 double corr = std::abs(*t - *temp); // for boundary with constant temperature this will be zero anyway
485 if (corr > err) err = corr;
486 if (*temp > maxT) maxT = *temp;
487 }
488 if (err > toterr) toterr = err;
489
490 ++loopno;
491 ++loop;
492
493 // show max correction
494 this->writelog(LOG_RESULT, "Loop {:d}({:d}): max(T) = {:.3f} K, error = {:g} K", loop, loopno, maxT, err);
495
496 } while ((!this->iter_params.converged || err > maxerr) && (loops == 0 || loop < loops));
497
498 outTemperature.fireChanged();
499 outHeatFlux.fireChanged();
500
501 return toterr;
502}
503
504 template<typename Geometry2DType>
505 void ThermalFem2DSolver<Geometry2DType>::saveHeatFluxes()
506{
507 this->writelog(LOG_DETAIL, "Computing heat fluxes");
508
509 fluxes.reset(this->maskedMesh->getElementsCount());
510
511 for (auto e: this->maskedMesh->elements())
512 {
513 Vec<2,double> midpoint = e.getMidpoint();
514 auto material = this->geometry->getMaterial(midpoint);
515
516 size_t loleftno = e.getLoLoIndex();
517 size_t lorghtno = e.getUpLoIndex();
518 size_t upleftno = e.getLoUpIndex();
519 size_t uprghtno = e.getUpUpIndex();
520
521 double temp = 0.25 * (temperatures[loleftno] + temperatures[lorghtno] +
522 temperatures[upleftno] + temperatures[uprghtno]);
523
524 double kx, ky;
525 auto leaf = dynamic_pointer_cast<const GeometryObjectD<2>>(
526 this->geometry->getMatchingAt(midpoint, &GeometryObject::PredicateIsLeaf)
527 );
528 if (leaf)
529 std::tie(kx,ky) = std::tuple<double,double>(material->thermk(temp, leaf->getBoundingBox().height()));
530 else
531 std::tie(kx,ky) = std::tuple<double,double>(material->thermk(temp));
532
533
534 fluxes[e.getIndex()] = vec(
535 - 0.5e6 * kx * (- temperatures[loleftno] + temperatures[lorghtno]
536 - temperatures[upleftno] + temperatures[uprghtno]) / (e.getUpper0() - e.getLower0()), // 1e6 - from um to m
537 - 0.5e6 * ky * (- temperatures[loleftno] - temperatures[lorghtno]
538 + temperatures[upleftno] + temperatures[uprghtno]) / (e.getUpper1() - e.getLower1())); // 1e6 - from um to m
539 }
540}
541
542
543 template<typename Geometry2DType>
544 const LazyData<double> ThermalFem2DSolver<Geometry2DType>::getTemperatures(const shared_ptr<const MeshD<2>>& dest_mesh, InterpolationMethod method) const {
545 this->writelog(LOG_DEBUG, "Getting temperatures");
546 if (!temperatures) return LazyData<double>(dest_mesh->size(), inittemp); // in case the receiver is connected and no temperature calculated yet
547 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
548 if (this->maskedMesh->full())
549 return SafeData<double>(interpolate(this->mesh, temperatures, dest_mesh, method, this->geometry), 300.);
550 else
551 return SafeData<double>(interpolate(this->maskedMesh, temperatures, dest_mesh, method, this->geometry), 300.);
552}
553
554
555 template<typename Geometry2DType>
556 const LazyData<Vec<2>> ThermalFem2DSolver<Geometry2DType>::getHeatFluxes(const shared_ptr<const MeshD<2>>& dest_mesh, InterpolationMethod method) {
557 this->writelog(LOG_DEBUG, "Getting heat fluxes");
558 if (!temperatures) return LazyData<Vec<2>>(dest_mesh->size(), Vec<2>(0.,0.)); // in case the receiver is connected and no fluxes calculated yet
559 if (!fluxes) saveHeatFluxes(); // we will compute fluxes only if they are needed
560 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
561 if (this->maskedMesh->full())
562 return SafeData<Vec<2>>(interpolate(this->mesh->getElementMesh(), fluxes, dest_mesh, method,
563 InterpolationFlags(this->geometry, InterpolationFlags::Symmetry::NP, InterpolationFlags::Symmetry::PN)),
564 Zero<Vec<2>>());
565 else
566 return SafeData<Vec<2>>(interpolate(this->maskedMesh->getElementMesh(), fluxes, dest_mesh, method,
567 InterpolationFlags(this->geometry, InterpolationFlags::Symmetry::NP, InterpolationFlags::Symmetry::PN)),
568 Zero<Vec<2>>());
569}
570
571
572 template<typename Geometry2DType> ThermalFem2DSolver<Geometry2DType>::
573 ThermalConductivityData::ThermalConductivityData(const ThermalFem2DSolver<Geometry2DType>* solver, const shared_ptr<const MeshD<2>>& dst_mesh):
574 solver(solver), dest_mesh(dst_mesh), flags(solver->geometry)
575{
576 if (solver->temperatures) temps = interpolate(solver->maskedMesh, solver->temperatures, solver->maskedMesh->getElementMesh(), INTERPOLATION_LINEAR);
577 else temps = LazyData<double>(solver->maskedMesh->getElementsCount(), solver->inittemp);
578}
579
580 template<typename Geometry2DType> Tensor2<double> ThermalFem2DSolver<Geometry2DType>::
581 ThermalConductivityData::at(std::size_t i) const {
582 auto point = flags.wrap(dest_mesh->at(i));
583 std::size_t x = solver->mesh->axis[0]->findUpIndex(point[0]),
584 y = solver->mesh->axis[1]->findUpIndex(point[1]);
585 if (x == 0 || y == 0 || x == solver->mesh->axis[0]->size() || y == solver->mesh->axis[1]->size())
586 return Tensor2<double>(NAN);
587 else {
588 auto elem = solver->maskedMesh->element(x-1, y-1);
589 size_t idx = elem.getIndex();
590 if (idx == RectangularMaskedMesh2D::Element::UNKNOWN_ELEMENT_INDEX) return Tensor2<double>(NAN);
591 auto material = solver->geometry->getMaterial(elem.getMidpoint());
592 return material->thermk(temps[idx], solver->thickness[idx]);
593 }
594}
595
596 template<typename Geometry2DType> std::size_t ThermalFem2DSolver<Geometry2DType>::
597 ThermalConductivityData::size() const { return dest_mesh->size(); }
598
599 template<typename Geometry2DType>
600 const LazyData<Tensor2<double>> ThermalFem2DSolver<Geometry2DType>::getThermalConductivity(const shared_ptr<const MeshD<2>>& dst_mesh, InterpolationMethod) {
601 this->initCalculation();
602 this->writelog(LOG_DEBUG, "Getting thermal conductivities");
603 return LazyData<Tensor2<double>>(new
604 ThermalFem2DSolver<Geometry2DType>::ThermalConductivityData(this, dst_mesh)
605 );
606}
607
608
609 template<> std::string ThermalFem2DSolver<Geometry2DCartesian>::getClassName() const { return "thermal.Static2D"; }
610 template<> std::string ThermalFem2DSolver<Geometry2DCylindrical>::getClassName() const { return "thermal.StaticCyl"; }
611
612 template struct PLASK_SOLVER_API ThermalFem2DSolver<Geometry2DCartesian>;
613 template struct PLASK_SOLVER_API ThermalFem2DSolver<Geometry2DCylindrical>;
614
615}}} // namespace plask::thermal::tstatic