PLaSK library
Loading...
Searching...
No Matches
femT2d.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 "femT2d.hpp"
15
16namespace plask { namespace thermal { namespace dynamic {
17
18const double BIG = 1e16;
19
20template<typename Geometry2DType>
21DynamicThermalFem2DSolver<Geometry2DType>::DynamicThermalFem2DSolver(const std::string& name) :
22 FemSolverWithMaskedMesh<Geometry2DType, RectangularMesh<2>>(name),
23 outTemperature(this, &DynamicThermalFem2DSolver<Geometry2DType>::getTemperatures),
24 outHeatFlux(this, &DynamicThermalFem2DSolver<Geometry2DType>::getHeatFluxes),
25 outThermalConductivity(this, &DynamicThermalFem2DSolver<Geometry2DType>::getThermalConductivity),
26 inittemp(300.),
27 methodparam(0.5),
28 timestep(0.1),
29 elapstime(0.),
30 lumping(true),
31 rebuildfreq(0),
32 logfreq(500)
33{
34 temperatures.reset();
35 fluxes.reset();
36 inHeat = 0.;
37}
38
39
40template<typename Geometry2DType>
43
44
45template<typename Geometry2DType>
46void DynamicThermalFem2DSolver<Geometry2DType>::loadConfiguration(XMLReader &source, Manager &manager)
47{
48 while (source.requireTagOrEnd())
49 {
50 std::string param = source.getNodeName();
51
52 if (param == "temperature")
53 this->readBoundaryConditions(manager, source, temperature_boundary);
54
55 else if (param == "loop") {
56 inittemp = source.getAttribute<double>("inittemp", inittemp);
57 timestep = source.getAttribute<double>("timestep", timestep);
58 rebuildfreq = source.getAttribute<size_t>("rebuildfreq", rebuildfreq);
59 logfreq = source.getAttribute<size_t>("logfreq", logfreq);
60 source.requireTagEnd();
61 }
62
63 else if (source.getNodeName() == "matrix") {
64 methodparam = source.getAttribute<double>("methodparam", methodparam);
65 lumping = source.getAttribute<bool>("lumping", lumping);
66 this->parseFemConfiguration(source, manager);
67 } else {
68 this->parseStandardConfiguration(source, manager);
69 }
70 }
71}
72
73
74template<typename Geometry2DType>
75void DynamicThermalFem2DSolver<Geometry2DType>::onInitialize() {
76 if (!this->geometry) throw NoGeometryException(this->getId());
77 if (!this->mesh) throw NoMeshException(this->getId());
78 elapstime = 0.;
79
80 FemSolverWithMaskedMesh<Geometry2DType, RectangularMesh<2>>::onInitialize();
81
82 temperatures.reset(this->maskedMesh->size(), inittemp);
83
84 thickness.reset(this->maskedMesh->getElementsCount(), NAN);
85 // Set stiffness matrix and load vector
86 for (auto elem: this->maskedMesh->elements())
87 {
88 if (!isnan(thickness[elem.getIndex()])) continue;
89 auto material = this->geometry->getMaterial(elem.getMidpoint());
90 double top = elem.getUpper1(), bottom = elem.getLower1();
91 size_t row = elem.getIndex1();
92 size_t itop = row+1, ibottom = row;
93 size_t c = elem.getIndex0();
94 for (size_t r = row; r > 0; r--) {
95 auto e = this->mesh->element(c, r-1);
96 auto m = this->geometry->getMaterial(e.getMidpoint());
97 if (m == material) { //TODO ignore doping
98 bottom = e.getLower1();
99 ibottom = r-1;
100 }
101 else break;
102 }
103 for (size_t r = elem.getIndex1()+1; r < this->mesh->axis[1]->size()-1; r++) {
104 auto e = this->mesh->element(c, r);
105 auto m = this->geometry->getMaterial(e.getMidpoint());
106 if (m == material) { //TODO ignore doping
107 top = e.getUpper1();
108 itop = r+1;
109 }
110 else break;
111 }
112 double h = top - bottom;
113 for (size_t r = ibottom; r != itop; ++r) {
114 size_t idx = this->maskedMesh->element(c, r).getIndex();
115 if (idx != RectangularMaskedMesh2D::Element::UNKNOWN_ELEMENT_INDEX)
116 thickness[idx] = h;
117 }
118 }
119}
120
121
122template<typename Geometry2DType> void DynamicThermalFem2DSolver<Geometry2DType>::onInvalidate() {
123 temperatures.reset();
124 fluxes.reset();
125}
126
127
128template<>
129void DynamicThermalFem2DSolver<Geometry2DCartesian>::setMatrix(
130 FemMatrix& A, FemMatrix& B, DataVector<double>& F,
131 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,double>& btemperature)
132{
133 this->writelog(LOG_DETAIL, "Setting up matrix system ({})", A.describe());
134
135 auto heatdensities = inHeat(this->maskedMesh->getElementMesh());
136
137 // zero the matrices A, B and the load vector F
138 A.clear();
139 B.clear();
140 F.fill(0.);
141
142 // Set stiffness matrix and load vector
143 for (auto elem: this->maskedMesh->elements())
144 {
145 // nodes numbers for the current element
146 size_t loleftno = elem.getLoLoIndex();
147 size_t lorghtno = elem.getUpLoIndex();
148 size_t upleftno = elem.getLoUpIndex();
149 size_t uprghtno = elem.getUpUpIndex();
150
151 // element size
152 double elemwidth = elem.getUpper0() - elem.getLower0();
153 double elemheight = elem.getUpper1() - elem.getLower1();
154
155 // point and material in the middle of the element
156 Vec<2,double> midpoint = elem.getMidpoint();
157 auto material = this->geometry->getMaterial(midpoint);
158
159 // average temperature on the element
160 double temp = 0.25 * (temperatures[loleftno] + temperatures[lorghtno] + temperatures[upleftno] + temperatures[uprghtno]);
161
162 // thermal conductivity
163 double kx, ky;
164 std::tie(kx,ky) = std::tuple<double,double>(material->thermk(temp, thickness[elem.getIndex()]));
165
166 // element of heat capacity matrix
167 double c = material->cp(temp) * material->dens(temp) * 0.25 * 1E-12 * elemheight * elemwidth / timestep / 1E-9;
168
169 kx *= elemheight; kx /= elemwidth;
170 ky *= elemwidth; ky /= elemheight;
171
172 // load vector: heat densities
173 double f = 0.25e-12 * elemwidth * elemheight * heatdensities[elem.getIndex()]; // 1e-12 -> to transform µm² into m²
174
175 // set symmetric matrix components in thermal conductivity matrix
176 double k44, k33, k22, k11, k43, k21, k42, k31, k32, k41;
177
178 k44 = k33 = k22 = k11 = (kx + ky) / 3.;
179 k43 = k21 = (-2. * kx + ky) / 6.;
180 k42 = k31 = - (kx + ky) / 6.;
181 k32 = k41 = (kx - 2. * ky) / 6.;
182
183 //Whether lumping the mass matrices A, B?
184 if (lumping)
185 {
186 A(loleftno, loleftno) += methodparam*k11 + c;
187 A(lorghtno, lorghtno) += methodparam*k22 + c;
188 A(uprghtno, uprghtno) += methodparam*k33 + c;
189 A(upleftno, upleftno) += methodparam*k44 + c;
190
191 A(lorghtno, loleftno) += methodparam*k21;
192 A(uprghtno, loleftno) += methodparam*k31;
193 A(upleftno, loleftno) += methodparam*k41;
194 A(uprghtno, lorghtno) += methodparam*k32;
195 A(upleftno, lorghtno) += methodparam*k42;
196 A(upleftno, uprghtno) += methodparam*k43;
197
198 B(loleftno, loleftno) += -(1-methodparam)*k11 + c;
199 B(lorghtno, lorghtno) += -(1-methodparam)*k22 + c;
200 B(uprghtno, uprghtno) += -(1-methodparam)*k33 + c;
201 B(upleftno, upleftno) += -(1-methodparam)*k44 + c;
202
203 B(lorghtno, loleftno) += -(1-methodparam)*k21;
204 B(uprghtno, loleftno) += -(1-methodparam)*k31;
205 B(upleftno, loleftno) += -(1-methodparam)*k41;
206 B(uprghtno, lorghtno) += -(1-methodparam)*k32;
207 B(upleftno, lorghtno) += -(1-methodparam)*k42;
208 B(upleftno, uprghtno) += -(1-methodparam)*k43;
209 }
210 else
211 {
212 A(loleftno, loleftno) += methodparam*k11 + 4./9.*c;
213 A(lorghtno, lorghtno) += methodparam*k22 + 4./9.*c;
214 A(uprghtno, uprghtno) += methodparam*k33 + 4./9.*c;
215 A(upleftno, upleftno) += methodparam*k44 + 4./9.*c;
216
217 A(lorghtno, loleftno) += methodparam*k21 + 2./9.*c;
218 A(uprghtno, loleftno) += methodparam*k31 + 1./9.*c;
219 A(upleftno, loleftno) += methodparam*k41 + 2./9.*c;
220 A(uprghtno, lorghtno) += methodparam*k32 + 2./9.*c;
221 A(upleftno, lorghtno) += methodparam*k42 + 1./9.*c;
222 A(upleftno, uprghtno) += methodparam*k43 + 2./9.*c;
223
224 B(loleftno, loleftno) += -(1-methodparam)*k11 + 4./9.*c;
225 B(lorghtno, lorghtno) += -(1-methodparam)*k22 + 4./9.*c;
226 B(uprghtno, uprghtno) += -(1-methodparam)*k33 + 4./9.*c;
227 B(upleftno, upleftno) += -(1-methodparam)*k44 + 4./9.*c;
228
229 B(lorghtno, loleftno) += -(1-methodparam)*k21 + 2./9.*c;
230 B(uprghtno, loleftno) += -(1-methodparam)*k31 + 1./9.*c;
231 B(upleftno, loleftno) += -(1-methodparam)*k41 + 2./9.*c;
232 B(uprghtno, lorghtno) += -(1-methodparam)*k32 + 2./9.*c;
233 B(upleftno, lorghtno) += -(1-methodparam)*k42 + 1./9.*c;
234 B(upleftno, uprghtno) += -(1-methodparam)*k43 + 2./9.*c;
235 }
236 // set load vector
237 F[loleftno] += f;
238 F[lorghtno] += f;
239 F[uprghtno] += f;
240 F[upleftno] += f;
241 }
242
243 A.applyBC(btemperature, F);
244
245 // macierz A -> L L^T
246 A.factorize();
247
248#ifndef NDEBUG
249 double* aend = A.data + A.size;
250 for (double* pa = A.data; pa != aend; ++pa) {
251 if (isnan(*pa) || isinf(*pa))
252 throw ComputationError(this->getId(), "error in stiffness matrix at position {0}", pa-A.data);
253 }
254#endif
255
256}
257
258template<>
259void DynamicThermalFem2DSolver<Geometry2DCylindrical>::setMatrix(
260 FemMatrix& A, FemMatrix& B, DataVector<double>& F,
261 const BoundaryConditionsWithMesh<RectangularMesh<2>::Boundary,double>& btemperature)
262{
263 this->writelog(LOG_DETAIL, "Setting up matrix system ({})", A.describe());
264
265 auto heatdensities = inHeat(this->maskedMesh->getElementMesh());
266
267 // zero the matrices A, B and the load vector F
268 A.clear();
269 B.clear();
270 F.fill(0.);
271
272 // Set stiffness matrix and load vector
273 for (auto elem: this->maskedMesh->elements())
274 {
275 // nodes numbers for the current element
276 size_t loleftno = elem.getLoLoIndex();
277 size_t lorghtno = elem.getUpLoIndex();
278 size_t upleftno = elem.getLoUpIndex();
279 size_t uprghtno = elem.getUpUpIndex();
280
281 // element size
282 double elemwidth = elem.getUpper0() - elem.getLower0();
283 double elemheight = elem.getUpper1() - elem.getLower1();
284
285 // point and material in the middle of the element
286 Vec<2,double> midpoint = elem.getMidpoint();
287 double r = midpoint.rad_r();
288 auto material = this->geometry->getMaterial(midpoint);
289
290 // average temperature on the element
291 double temp = 0.25 * (temperatures[loleftno] + temperatures[lorghtno] + temperatures[upleftno] + temperatures[uprghtno]);
292
293 // thermal conductivity
294 double kx, ky;
295 std::tie(kx,ky) = std::tuple<double,double>(material->thermk(temp, thickness[elem.getIndex()]));
296
297 // element of heat capacity matrix
298 double c = material->cp(temp) * material->dens(temp) * 0.25 * 1E-12 * r * elemheight * elemwidth / timestep / 1E-9;
299
300 kx *= elemheight; kx /= elemwidth; kx *= r;
301 ky *= elemwidth; ky /= elemheight; ky *= r;
302
303 // load vector: heat densities
304 double f = 0.25e-12 * r * elemwidth * elemheight * heatdensities[elem.getIndex()]; // 1e-12 -> to transform µm² into m²
305
306 // set symmetric matrix components in thermal conductivity matrix
307 double k44, k33, k22, k11, k43, k21, k42, k31, k32, k41;
308
309 k44 = k33 = k22 = k11 = (kx + ky) / 3.;
310 k43 = k21 = (-2. * kx + ky) / 6.;
311 k42 = k31 = - (kx + ky) / 6.;
312 k32 = k41 = (kx - 2. * ky) / 6.;
313
314 //Whether lumping the mass matrices A, B?
315 if (lumping) {
316 A(loleftno, loleftno) += methodparam*k11 + c;
317 A(lorghtno, lorghtno) += methodparam*k22 + c;
318 A(uprghtno, uprghtno) += methodparam*k33 + c;
319 A(upleftno, upleftno) += methodparam*k44 + c;
320
321 A(lorghtno, loleftno) += methodparam*k21;
322 A(uprghtno, loleftno) += methodparam*k31;
323 A(upleftno, loleftno) += methodparam*k41;
324 A(uprghtno, lorghtno) += methodparam*k32;
325 A(upleftno, lorghtno) += methodparam*k42;
326 A(upleftno, uprghtno) += methodparam*k43;
327
328 B(loleftno, loleftno) += -(1-methodparam)*k11 + c;
329 B(lorghtno, lorghtno) += -(1-methodparam)*k22 + c;
330 B(uprghtno, uprghtno) += -(1-methodparam)*k33 + c;
331 B(upleftno, upleftno) += -(1-methodparam)*k44 + c;
332
333 B(lorghtno, loleftno) += -(1-methodparam)*k21;
334 B(uprghtno, loleftno) += -(1-methodparam)*k31;
335 B(upleftno, loleftno) += -(1-methodparam)*k41;
336 B(uprghtno, lorghtno) += -(1-methodparam)*k32;
337 B(upleftno, lorghtno) += -(1-methodparam)*k42;
338 B(upleftno, uprghtno) += -(1-methodparam)*k43;
339 } else {
340 A(loleftno, loleftno) += methodparam*k11 + 4./9.*c;
341 A(lorghtno, lorghtno) += methodparam*k22 + 4./9.*c;
342 A(uprghtno, uprghtno) += methodparam*k33 + 4./9.*c;
343 A(upleftno, upleftno) += methodparam*k44 + 4./9.*c;
344
345 A(lorghtno, loleftno) += methodparam*k21 + 2./9.*c;
346 A(uprghtno, loleftno) += methodparam*k31 + 1./9.*c;
347 A(upleftno, loleftno) += methodparam*k41 + 2./9.*c;
348 A(uprghtno, lorghtno) += methodparam*k32 + 2./9.*c;
349 A(upleftno, lorghtno) += methodparam*k42 + 1./9.*c;
350 A(upleftno, uprghtno) += methodparam*k43 + 2./9.*c;
351
352 B(loleftno, loleftno) += -(1-methodparam)*k11 + 4./9.*c;
353 B(lorghtno, lorghtno) += -(1-methodparam)*k22 + 4./9.*c;
354 B(uprghtno, uprghtno) += -(1-methodparam)*k33 + 4./9.*c;
355 B(upleftno, upleftno) += -(1-methodparam)*k44 + 4./9.*c;
356
357 B(lorghtno, loleftno) += -(1-methodparam)*k21 + 2./9.*c;
358 B(uprghtno, loleftno) += -(1-methodparam)*k31 + 1./9.*c;
359 B(upleftno, loleftno) += -(1-methodparam)*k41 + 2./9.*c;
360 B(uprghtno, lorghtno) += -(1-methodparam)*k32 + 2./9.*c;
361 B(upleftno, lorghtno) += -(1-methodparam)*k42 + 1./9.*c;
362 B(upleftno, uprghtno) += -(1-methodparam)*k43 + 2./9.*c;
363 }
364 // set load vector
365 F[loleftno] += f;
366 F[lorghtno] += f;
367 F[uprghtno] += f;
368 F[upleftno] += f;
369 }
370
371 A.applyBC(btemperature, F);
372
373 // macierz A -> L L^T
374 A.factorize();
375
376#ifndef NDEBUG
377 double* aend = A.data + A.size;
378 for (double* pa = A.data; pa != aend; ++pa) {
379 if (isnan(*pa) || isinf(*pa))
380 throw ComputationError(this->getId(), "error in stiffness matrix at position {0}", pa-A.data);
381 }
382#endif
383
384}
385
386
387template<typename Geometry2DType>
388double DynamicThermalFem2DSolver<Geometry2DType>::compute(double time)
389{
390 this->initCalculation();
391
392 fluxes.reset();
393
394 // store boundary conditions for current mesh
395 auto btemperature = temperature_boundary(this->maskedMesh, this->geometry);
396
397 size_t size = this->maskedMesh->size();
398
399 std::unique_ptr<FemMatrix> pA(this->getMatrix());
400 FemMatrix& A = *pA.get();
401 std::unique_ptr<FemMatrix> pB(this->getMatrix());
402 FemMatrix& B = *pB.get();
403
404 this->writelog(LOG_INFO, "Running thermal calculations");
405 maxT = *std::max_element(temperatures.begin(), temperatures.end());
406
407# ifndef NDEBUG
408 if (!temperatures.unique()) this->writelog(LOG_DEBUG, "Temperature data held by something else...");
409# endif
410 temperatures = temperatures.claim();
411 DataVector<double> F(size), X(size);
412
413 setMatrix(A, B, F, btemperature);
414
415 size_t r = rebuildfreq,
416 l = logfreq;
417
418 time += timestep/2.;
419 for (double t = 0.; t < time; t += timestep) {
420
421 if (rebuildfreq && r == 0)
422 {
423 setMatrix(A, B, F, btemperature);
424 r = rebuildfreq;
425 }
426
427 B.mult(temperatures, X);
428 for (std::size_t i = 0; i < X.size(); ++i) X[i] += F[i];
429
430 DataVector<double> T(X);
431
432 A.solverhs(T, temperatures);
433
434 if (logfreq && l == 0)
435 {
436 maxT = *std::max_element(temperatures.begin(), temperatures.end());
437 this->writelog(LOG_RESULT, "Time {:.2f} ns: max(T) = {:.3f} K", elapstime, maxT);
438 l = logfreq;
439 }
440
441 r--;
442 l--;
443 elapstime += timestep;
444 }
445
446 elapstime -= timestep;
447 outTemperature.fireChanged();
448 outHeatFlux.fireChanged();
449
450 return 0.;
451}
452
453template<typename Geometry2DType>
454void DynamicThermalFem2DSolver<Geometry2DType>::saveHeatFluxes()
455{
456 this->writelog(LOG_DETAIL, "Computing heat fluxes");
457
458 fluxes.reset(this->maskedMesh->getElementsCount());
459
460 for (auto e: this->maskedMesh->elements())
461 {
462 Vec<2,double> midpoint = e.getMidpoint();
463 auto material = this->geometry->getMaterial(midpoint);
464
465 size_t loleftno = e.getLoLoIndex();
466 size_t lorghtno = e.getUpLoIndex();
467 size_t upleftno = e.getLoUpIndex();
468 size_t uprghtno = e.getUpUpIndex();
469
470 double temp = 0.25 * (temperatures[loleftno] + temperatures[lorghtno] +
471 temperatures[upleftno] + temperatures[uprghtno]);
472
473 double kx, ky;
474 auto leaf = dynamic_pointer_cast<const GeometryObjectD<2>>(
475 this->geometry->getMatchingAt(midpoint, &GeometryObject::PredicateIsLeaf)
476 );
477 if (leaf)
478 std::tie(kx,ky) = std::tuple<double,double>(material->thermk(temp, leaf->getBoundingBox().height()));
479 else
480 std::tie(kx,ky) = std::tuple<double,double>(material->thermk(temp));
481
482
483 fluxes[e.getIndex()] = vec(
484 - 0.5e6 * kx * (- temperatures[loleftno] + temperatures[lorghtno]
485 - temperatures[upleftno] + temperatures[uprghtno]) / (e.getUpper0() - e.getLower0()), // 1e6 - from um to m
486 - 0.5e6 * ky * (- temperatures[loleftno] - temperatures[lorghtno]
487 + temperatures[upleftno] + temperatures[uprghtno]) / (e.getUpper1() - e.getLower1())); // 1e6 - from um to m
488 }
489}
490
491
492template<typename Geometry2DType>
493const LazyData<double> DynamicThermalFem2DSolver<Geometry2DType>::getTemperatures(const shared_ptr<const MeshD<2>>& dest_mesh, InterpolationMethod method) const {
494 this->writelog(LOG_DEBUG, "Getting temperatures");
495 if (!temperatures) return LazyData<double>(dest_mesh->size(), inittemp); // in case the receiver is connected and no temperature calculated yet
496 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
497 if (this->maskedMesh->full())
498 return SafeData<double>(interpolate(this->mesh, temperatures, dest_mesh, method, this->geometry), 300.);
499 else
500 return SafeData<double>(interpolate(this->maskedMesh, temperatures, dest_mesh, method, this->geometry), 300.);
501}
502
503
504template<typename Geometry2DType>
505const LazyData<Vec<2>> DynamicThermalFem2DSolver<Geometry2DType>::getHeatFluxes(const shared_ptr<const MeshD<2>>& dest_mesh, InterpolationMethod method) {
506 this->writelog(LOG_DEBUG, "Getting heat fluxes");
507 if (!temperatures) return LazyData<Vec<2>>(dest_mesh->size(), Vec<2>(0.,0.)); // in case the receiver is connected and no fluxes calculated yet
508 if (!fluxes) saveHeatFluxes(); // we will compute fluxes only if they are needed
509 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
510 if (this->maskedMesh->full())
511 return SafeData<Vec<2>>(interpolate(this->mesh->getElementMesh(), fluxes, dest_mesh, method,
512 InterpolationFlags(this->geometry, InterpolationFlags::Symmetry::NP, InterpolationFlags::Symmetry::PN)),
513 Zero<Vec<2>>());
514 else
515 return SafeData<Vec<2>>(interpolate(this->maskedMesh->getElementMesh(), fluxes, dest_mesh, method,
516 InterpolationFlags(this->geometry, InterpolationFlags::Symmetry::NP, InterpolationFlags::Symmetry::PN)),
517 Zero<Vec<2>>());
518}
519
520
521template<typename Geometry2DType> DynamicThermalFem2DSolver<Geometry2DType>::
522ThermalConductivityData::ThermalConductivityData(const DynamicThermalFem2DSolver<Geometry2DType>* solver, const shared_ptr<const MeshD<2>>& dst_mesh):
523 solver(solver), dest_mesh(dst_mesh), flags(solver->geometry)
524{
525 if (solver->temperatures) temps = interpolate(solver->maskedMesh, solver->temperatures, solver->maskedMesh->getElementMesh(), INTERPOLATION_LINEAR);
526 else temps = LazyData<double>(solver->maskedMesh->getElementsCount(), solver->inittemp);
527}
528
529template<typename Geometry2DType> Tensor2<double> DynamicThermalFem2DSolver<Geometry2DType>::
530ThermalConductivityData::at(std::size_t i) const {
531 auto point = flags.wrap(dest_mesh->at(i));
532 size_t x = solver->mesh->axis[0]->findUpIndex(point[0]),
533 y = solver->mesh->axis[1]->findUpIndex(point[1]);
534 if (x == 0 || y == 0 || x == solver->mesh->axis[0]->size() || y == solver->mesh->axis[1]->size())
535 return Tensor2<double>(NAN);
536 else {
537 auto elem = solver->maskedMesh->element(x-1, y-1);
538 auto material = solver->geometry->getMaterial(elem.getMidpoint());
539 size_t idx = elem.getIndex();
540 if (idx == RectangularMaskedMesh2D::Element::UNKNOWN_ELEMENT_INDEX) return Tensor2<double>(NAN);
541 return material->thermk(temps[idx], solver->thickness[idx]);
542 }
543}
544
545template<typename Geometry2DType> std::size_t DynamicThermalFem2DSolver<Geometry2DType>::
546ThermalConductivityData::size() const { return dest_mesh->size(); }
547
548template<typename Geometry2DType>
556
557
558template<> std::string DynamicThermalFem2DSolver<Geometry2DCartesian>::getClassName() const { return "thermal.Dynamic2D"; }
559template<> std::string DynamicThermalFem2DSolver<Geometry2DCylindrical>::getClassName() const { return "thermal.DynamicCyl"; }
560
561template struct PLASK_SOLVER_API DynamicThermalFem2DSolver<Geometry2DCartesian>;
562template struct PLASK_SOLVER_API DynamicThermalFem2DSolver<Geometry2DCylindrical>;
563
564}}} // namespace plask::thermal::thermal