PLaSK library

Loading...
Searching...
No Matches
electr3d.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 <type_traits>
15
16#include "electr3d.hpp"
17
18 namespace plask { namespace electrical { namespace shockley {
19
20 ElectricalFem3DSolver::ElectricalFem3DSolver(const std::string& name)
21 : FemSolverWithMaskedMesh<Geometry3D, plask::RectangularMesh<3>>(name),
22 pcond(5.),
23 ncond(50.),
24 loopno(0),
25 default_junction_conductivity(Tensor2<double>(0., 5.)),
26 maxerr(0.05),
27 outVoltage(this, &ElectricalFem3DSolver::getVoltage),
28 outCurrentDensity(this, &ElectricalFem3DSolver::getCurrentDensity),
29 outHeat(this, &ElectricalFem3DSolver::getHeatDensity),
30 outConductivity(this, &ElectricalFem3DSolver::getConductivity),
31 convergence(CONVERGENCE_FAST) {
32 potential.reset();
33 current.reset();
34 inTemperature = 300.;
35 junction_conductivity.reset(1, default_junction_conductivity);
36 algorithm = ALGORITHM_ITERATIVE;
37}
38
39 ElectricalFem3DSolver::~ElectricalFem3DSolver() {}
40
41 void ElectricalFem3DSolver::loadConfiguration(XMLReader& source, Manager& manager) {
42 while (source.requireTagOrEnd()) parseConfiguration(source, manager);
43}
44
45 void ElectricalFem3DSolver::parseConfiguration(XMLReader& source, Manager& manager) {
46 std::string param = source.getNodeName();
47
48 if (param == "voltage")
49 readBoundaryConditions(manager, source, voltage_boundary);
50
51 else if (param == "loop") {
52 if (source.hasAttribute("start-cond") || source.hasAttribute("start-cond-inplane")) {
53 double c0 = default_junction_conductivity.c00, c1 = default_junction_conductivity.c11;
54 auto vert = source.getAttribute<std::string>("start-cond");
55 if (vert) {
56 if (vert->find(',') == std::string::npos) {
57 c1 = boost::lexical_cast<double>(*vert);
58 } else {
59 if (source.hasAttribute("start-cond-inplane"))
60 throw XMLException(
61 source, "tag attribute 'start-cond' has two values, but attribute 'start-cond-inplane' is also provided");
62 auto values = splitString2(*vert, ',');
63 c0 = boost::lexical_cast<double>(values.first);
64 c1 = boost::lexical_cast<double>(values.second);
65 }
66 }
67 c0 = source.getAttribute<double>("start-cond-inplane", c0);
68 this->setCondJunc(Tensor2<double>(c0, c1));
69 }
70 convergence = source.enumAttribute<Convergence>("convergence")
71 .value("fast", CONVERGENCE_FAST)
72 .value("stable", CONVERGENCE_STABLE)
73 .get(convergence);
74 maxerr = source.getAttribute<double>("maxerr", maxerr);
75 source.requireTagEnd();
76 }
77
78 else if (param == "contacts") {
79 pcond = source.getAttribute<double>("pcond", pcond);
80 ncond = source.getAttribute<double>("ncond", ncond);
81 source.requireTagEnd();
82 }
83
84 else if (!this->parseFemConfiguration(source, manager)) {
85 this->parseStandardConfiguration(source, manager);
86 }
87}
88
89 void ElectricalFem3DSolver::setupActiveRegions() {
90 if (!geometry || !mesh) {
91 if (junction_conductivity.size() != 1) {
92 Tensor2<double> condy(0., 0.);
93 for (auto cond : junction_conductivity) condy += cond;
94 junction_conductivity.reset(1, condy / double(junction_conductivity.size()));
95 }
96 return;
97 }
98
99 setupMaskedMesh();
100
101 shared_ptr<RectangularMesh<3>> points = mesh->getElementMesh();
102
103 std::map<size_t, Active::Region> regions;
104 size_t nreg = 0;
105
106 for (size_t lon = 0; lon < points->axis[0]->size(); ++lon) {
107 for (size_t tra = 0; tra < points->axis[1]->size(); ++tra) {
108 size_t num = 0;
109 size_t start = 0;
110 for (size_t ver = 0; ver < points->axis[2]->size(); ++ver) {
111 auto point = points->at(lon, tra, ver);
112 size_t cur = isActive(point);
113 if (cur != num) {
114 if (num) { // summarize current region
115 auto found = regions.find(num);
116 if (found == regions.end()) { // `num` is a new region
117 regions[num] = Active::Region(start, ver, lon, tra);
118 if (nreg < num) nreg = num;
119 } else {
120 Active::Region& region = found->second;
121 if (start != region.bottom || ver != region.top)
122 throw Exception("{0}: Junction {1} does not have top and bottom edges at constant heights",
123 this->getId(), num - 1);
124 if (tra < region.left) region.left = tra;
125 if (tra >= region.right) region.right = tra + 1;
126 if (lon < region.back) region.back = lon;
127 if (lon >= region.front) region.front = lon + 1;
128 }
129 }
130 num = cur;
131 start = ver;
132 }
133 if (cur) {
134 auto found = regions.find(cur);
135 if (found != regions.end()) {
136 Active::Region& region = found->second;
137 if (region.warn && lon != region.lon && tra != region.tra &&
138 *this->geometry->getMaterial(points->at(lon, tra, ver)) !=
139 *this->geometry->getMaterial(points->at(region.lon, region.tra, ver))) {
140 writelog(LOG_WARNING, "Junction {} is laterally non-uniform", num - 1);
141 region.warn = false;
142 }
143 }
144 }
145 }
146 if (num) { // summarize current region
147 auto found = regions.find(num);
148 if (found == regions.end()) { // `current` is a new region
149 regions[num] = Active::Region(start, points->axis[2]->size(), lon, tra);
150 } else {
151 Active::Region& region = found->second;
152 if (start != region.bottom || points->axis[2]->size() != region.top)
153 throw Exception("{0}: Junction {1} does not have top and bottom edges at constant heights", this->getId(),
154 num - 1);
155 if (tra < region.left) region.left = tra;
156 if (tra >= region.right) region.right = tra + 1;
157 if (lon < region.back) region.back = lon;
158 if (lon >= region.front) region.front = lon + 1;
159 }
160 }
161 }
162 }
163
164 size_t condsize = 0;
165 active.resize(nreg);
166
167 for (auto& ireg : regions) {
168 size_t num = ireg.first - 1;
169 Active::Region& reg = ireg.second;
170 double height = this->mesh->axis[2]->at(reg.top) - this->mesh->axis[2]->at(reg.bottom);
171 active[num] = Active(condsize, reg, height);
172 condsize += (reg.right - reg.left) * (reg.front - reg.back);
173 this->writelog(LOG_DETAIL, "Detected junction {0} thickness = {1}nm", num, 1e3 * height);
174 this->writelog(LOG_DEBUG, "Junction {0} span: [{1},{3},{5}]-[{2},{4},{6}]", num, reg.back, reg.front, reg.left, reg.right,
175 reg.bottom, reg.top);
176 }
177
178 if (junction_conductivity.size() != condsize) {
179 Tensor2<double> condy(0., 0.);
180 for (auto cond : junction_conductivity) condy += cond;
181 junction_conductivity.reset(condsize, condy / double(junction_conductivity.size()));
182 }
183}
184
185 void ElectricalFem3DSolver::onInitialize() {
186 if (!geometry) throw NoGeometryException(getId());
187 if (!mesh) throw NoMeshException(getId());
188 setupActiveRegions();
189 loopno = 0;
190 potential.reset(maskedMesh->size(), 0.);
191 current.reset(maskedMesh->getElementsCount(), vec(0., 0., 0.));
192 conds.reset(maskedMesh->getElementsCount());
193}
194
195 void ElectricalFem3DSolver::onInvalidate() {
196 conds.reset();
197 potential.reset();
198 current.reset();
199 heat.reset();
200 junction_conductivity.reset(1, default_junction_conductivity);
201}
202
203 LazyData<double> ElectricalFem3DSolver::loadConductivity() {
204 auto midmesh = (this->maskedMesh)->getElementMesh();
205 auto temperature = inTemperature(midmesh);
206
207 for (auto e : this->maskedMesh->elements()) {
208 size_t i = e.getIndex();
209 Vec<3, double> midpoint = e.getMidpoint();
210
211 auto roles = this->geometry->getRolesAt(midpoint);
212 if (size_t actn = isActive(midpoint)) {
213 const auto& act = active[actn - 1];
214 conds[i] = junction_conductivity[act.offset + act.ld * e.getIndex1() + e.getIndex0()];
215 if (isnan(conds[i].c11) || abs(conds[i].c11) < 1e-16) conds[i].c11 = 1e-16;
216 } else if (roles.find("p-contact") != roles.end()) {
217 conds[i] = Tensor2<double>(pcond, pcond);
218 } else if (roles.find("n-contact") != roles.end()) {
219 conds[i] = Tensor2<double>(ncond, ncond);
220 } else
221 conds[i] = this->geometry->getMaterial(midpoint)->cond(temperature[i]);
222 }
223
224 return temperature;
225}
226
227 void ElectricalFem3DSolver::saveConductivity() {
228 for (size_t n = 0; n < active.size(); ++n) {
229 const auto& act = active[n];
230 size_t v = (act.top + act.bottom) / 2;
231 for (size_t t = act.left; t != act.right; ++t) {
232 size_t offset = act.offset + act.ld * t;
233 for (size_t l = act.back; l != act.front; ++l)
234 junction_conductivity[offset + l] = conds[this->maskedMesh->element(l, t, v).getIndex()];
235 }
236 }
237}
238
239 void ElectricalFem3DSolver::setMatrix(FemMatrix& A,
240 DataVector<double>& B,
241 const BoundaryConditionsWithMesh<RectangularMesh<3>::Boundary, double>& bvoltage,
242 const LazyData<double>& temperature) {
243 this->writelog(LOG_DETAIL, "Setting up matrix system ({})", A.describe());
244
245 // Update junction conductivities
246 if (loopno != 0) {
247 for (auto elem : maskedMesh->elements()) {
248 if (size_t nact = isActive(elem)) {
249 size_t index = elem.getIndex(), lll = elem.getLoLoLoIndex(), uuu = elem.getUpUpUpIndex();
250 size_t back = maskedMesh->index0(lll), front = maskedMesh->index0(uuu), left = maskedMesh->index1(lll),
251 right = maskedMesh->index1(uuu);
252 const Active& act = active[nact - 1];
253 double U =
254 0.25 *
255 (-potential[maskedMesh->index(back, left, act.bottom)] - potential[maskedMesh->index(front, left, act.bottom)] -
256 potential[maskedMesh->index(back, right, act.bottom)] -
257 potential[maskedMesh->index(front, right, act.bottom)] + potential[maskedMesh->index(back, left, act.top)] +
258 potential[maskedMesh->index(front, left, act.top)] + potential[maskedMesh->index(back, right, act.top)] +
259 potential[maskedMesh->index(front, right, act.top)]);
260 double jy = 0.1 * conds[index].c11 * U / act.height; // [j] = kA/cm²
261 size_t tidx = this->maskedMesh->element(elem.getIndex0(), elem.getIndex1(), (act.top + act.bottom) / 2).getIndex();
262 Tensor2<double> cond = activeCond(nact - 1, U, jy, temperature[tidx]);
263 switch (convergence) {
264 case CONVERGENCE_STABLE:
265 cond = 0.5 * (conds[index] + cond);
266 case CONVERGENCE_FAST:
267 conds[index] = cond;
268 }
269 if (isnan(conds[index].c11) || abs(conds[index].c11) < 1e-16) {
270 conds[index].c11 = 1e-16;
271 }
272 }
273 }
274 }
275
276 // Zero the matrix and the load vector
277 A.clear();
278 B.fill(0.);
279
280 // Set stiffness matrix and load vector
281 for (auto elem : maskedMesh->elements()) {
282 size_t index = elem.getIndex();
283
284 // nodes numbers for the current element
285 size_t idx[8];
286 idx[0] = elem.getLoLoLoIndex(); // z 4-----6
287 idx[1] = elem.getUpLoLoIndex(); // |__y /| /|
288 idx[2] = elem.getLoUpLoIndex(); // x/ 5-----7 |
289 idx[3] = elem.getUpUpLoIndex(); // | 0---|-2
290 idx[4] = elem.getLoLoUpIndex(); // |/ |/
291 idx[5] = elem.getUpLoUpIndex(); // 1-----3
292 idx[6] = elem.getLoUpUpIndex(); //
293 idx[7] = elem.getUpUpUpIndex(); //
294
295 // element size
296 double dx = elem.getUpper0() - elem.getLower0();
297 double dy = elem.getUpper1() - elem.getLower1();
298 double dz = elem.getUpper2() - elem.getLower2();
299
300 // point and material in the middle of the element
301 Vec<3> middle = elem.getMidpoint();
302 auto material = geometry->getMaterial(middle);
303
304 // average voltage on the element
305 double temp = 0.;
306 for (int i = 0; i < 8; ++i) temp += potential[idx[i]];
307 temp *= 0.125;
308
309 // electrical conductivity
310 double kx, ky = conds[index].c00, kz = conds[index].c11;
311
312 ky *= 1e-6;
313 kz *= 1e-6; // 1/m -> 1/µm
314 kx = ky;
315
316 kx /= dx;
317 kx *= dy;
318 kx *= dz;
319 ky *= dx;
320 ky /= dy;
321 ky *= dz;
322 kz *= dx;
323 kz *= dy;
324 kz /= dz;
325
326 // set symmetric matrix components
327 double K[8][8];
328 K[0][0] = K[1][1] = K[2][2] = K[3][3] = K[4][4] = K[5][5] = K[6][6] = K[7][7] = (kx + ky + kz) / 9.;
329
330 K[1][0] = K[3][2] = K[5][4] = K[7][6] = (-2. * kx + ky + kz) / 18.;
331 K[2][0] = K[3][1] = K[6][4] = K[7][5] = (kx - 2. * ky + kz) / 18.;
332 K[4][0] = K[5][1] = K[6][2] = K[7][3] = (kx + ky - 2. * kz) / 18.;
333
334 K[4][2] = K[5][3] = K[6][0] = K[7][1] = (kx - 2. * ky - 2. * kz) / 36.;
335 K[4][1] = K[5][0] = K[6][3] = K[7][2] = (-2. * kx + ky - 2. * kz) / 36.;
336 K[2][1] = K[3][0] = K[6][5] = K[7][4] = (-2. * kx - 2. * ky + kz) / 36.;
337
338 K[4][3] = K[5][2] = K[6][1] = K[7][0] = -(kx + ky + kz) / 36.;
339
340 for (int i = 0; i < 8; ++i)
341 for (int j = 0; j <= i; ++j) A(idx[i], idx[j]) += K[i][j];
342 }
343
344 A.applyBC(bvoltage, B);
345
346#ifndef NDEBUG
347 double* aend = A.data + A.size;
348 for (double* pa = A.data; pa != aend; ++pa) {
349 if (isnan(*pa) || isinf(*pa))
350 throw ComputationError(getId(), "error in stiffness matrix at position {0} ({1})", pa - A.data,
351 isnan(*pa) ? "nan" : "inf");
352 }
353#endif
354}
355
356 double ElectricalFem3DSolver::compute(unsigned loops) {
357 this->initCalculation();
358
359 // store boundary conditions for current mesh
360 auto bvoltage = voltage_boundary(maskedMesh, geometry);
361
362 this->writelog(LOG_INFO, "Running electrical calculations");
363
364 unsigned loop = 0;
365 double err = 0.;
366 toterr = 0.;
367
368 std::unique_ptr<FemMatrix> pA(this->getMatrix());
369 FemMatrix& A = *pA.get();
370
371#ifndef NDEBUG
372 if (!potential.unique()) this->writelog(LOG_DEBUG, "Potentials data held by something else...");
373#endif
374 potential = potential.claim();
375
376 DataVector<double> rhs(potential.size());
377
378 auto temperature = loadConductivity();
379
380 bool noactive = (active.size() == 0);
381 double minj = 100e-7; // assume no significant heating below this current
382
383 do {
384 setMatrix(A, rhs, bvoltage, temperature);
385 A.solve(rhs, potential);
386
387 err = 0.;
388 double mcur = 0.;
389 for (auto el : maskedMesh->elements()) {
390 size_t i = el.getIndex();
391 size_t lll = el.getLoLoLoIndex();
392 size_t llu = el.getLoLoUpIndex();
393 size_t lul = el.getLoUpLoIndex();
394 size_t luu = el.getLoUpUpIndex();
395 size_t ull = el.getUpLoLoIndex();
396 size_t ulu = el.getUpLoUpIndex();
397 size_t uul = el.getUpUpLoIndex();
398 size_t uuu = el.getUpUpUpIndex();
399 auto cur = vec(-0.025 * conds[i].c00 *
400 (-potential[lll] - potential[llu] - potential[lul] - potential[luu] + potential[ull] +
401 potential[ulu] + potential[uul] + potential[uuu]) /
402 (el.getUpper0() - el.getLower0()), // [j] = kA/cm²
403 -0.025 * conds[i].c00 *
404 (-potential[lll] - potential[llu] + potential[lul] + potential[luu] - potential[ull] -
405 potential[ulu] + potential[uul] + potential[uuu]) /
406 (el.getUpper1() - el.getLower1()), // [j] = kA/cm²
407 -0.025 * conds[i].c11 *
408 (-potential[lll] + potential[llu] - potential[lul] + potential[luu] - potential[ull] +
409 potential[ulu] - potential[uul] + potential[uuu]) /
410 (el.getUpper2() - el.getLower2()) // [j] = kA/cm²
411 );
412 if (noactive || isActive(el)) {
413 double acur = abs2(cur);
414 if (acur > mcur) {
415 mcur = acur;
416 maxcur = cur;
417 }
418 }
419 double delta = abs2(current[i] - cur);
420 if (delta > err) err = delta;
421 current[i] = cur;
422 }
423 mcur = sqrt(mcur);
424 err = 100. * sqrt(err) / max(mcur, minj);
425 if ((loop != 0 || mcur >= minj) && err > toterr) toterr = err;
426
427 ++loopno;
428 ++loop;
429
430 this->writelog(LOG_RESULT, "Loop {:d}({:d}): max(j{}) = {:g} kA/cm2, error = {:g}%", loop, loopno, noactive ? "" : "@junc",
431 mcur, err);
432
433 } while ((!iter_params.converged || err > maxerr) && (loops == 0 || loop < loops));
434
435 saveConductivity();
436
437 outVoltage.fireChanged();
438 outCurrentDensity.fireChanged();
439 outHeat.fireChanged();
440
441 return toterr;
442}
443
444 void ElectricalFem3DSolver::saveHeatDensity() {
445 this->writelog(LOG_DETAIL, "Computing heat densities");
446
447 heat.reset(maskedMesh->getElementsCount());
448
449 for (auto el : maskedMesh->elements()) {
450 size_t i = el.getIndex();
451 size_t lll = el.getLoLoLoIndex();
452 size_t llu = el.getLoLoUpIndex();
453 size_t lul = el.getLoUpLoIndex();
454 size_t luu = el.getLoUpUpIndex();
455 size_t ull = el.getUpLoLoIndex();
456 size_t ulu = el.getUpLoUpIndex();
457 size_t uul = el.getUpUpLoIndex();
458 size_t uuu = el.getUpUpUpIndex();
459 double dvx = -0.25e6 *
460 (-potential[lll] - potential[llu] - potential[lul] - potential[luu] + potential[ull] + potential[ulu] +
461 potential[uul] + potential[uuu]) /
462 (el.getUpper0() - el.getLower0()); // 1e6 - from µm to m
463 double dvy = -0.25e6 *
464 (-potential[lll] - potential[llu] + potential[lul] + potential[luu] - potential[ull] - potential[ulu] +
465 potential[uul] + potential[uuu]) /
466 (el.getUpper1() - el.getLower1()); // 1e6 - from µm to m
467 double dvz = -0.25e6 *
468 (-potential[lll] + potential[llu] - potential[lul] + potential[luu] - potential[ull] + potential[ulu] -
469 potential[uul] + potential[uuu]) /
470 (el.getUpper2() - el.getLower2()); // 1e6 - from µm to m
471 auto midpoint = el.getMidpoint();
472 if (geometry->getMaterial(midpoint)->kind() == Material::EMPTY || geometry->hasRoleAt("noheat", midpoint))
473 heat[i] = 0.;
474 else {
475 heat[i] = conds[i].c00 * dvx * dvx + conds[i].c00 * dvy * dvy + conds[i].c11 * dvz * dvz;
476 }
477 }
478}
479
480 double ElectricalFem3DSolver::integrateCurrent(size_t vindex, bool onlyactive) {
481 if (!potential) throw NoValue("current densities");
482 this->writelog(LOG_DETAIL, "Computing total current");
483 double result = 0.;
484 for (size_t i = 0; i < mesh->axis[0]->size() - 1; ++i) {
485 for (size_t j = 0; j < mesh->axis[1]->size() - 1; ++j) {
486 auto element = maskedMesh->element(i, j, vindex);
487 if (!onlyactive || isActive(element.getMidpoint())) {
488 size_t index = element.getIndex();
489 if (index != RectangularMaskedMesh3D::Element::UNKNOWN_ELEMENT_INDEX)
490 result += current[index].c2 * element.getSize0() * element.getSize1();
491 }
492 }
493 }
494 if (geometry->isSymmetric(Geometry::DIRECTION_LONG)) result *= 2.;
495 if (geometry->isSymmetric(Geometry::DIRECTION_TRAN)) result *= 2.;
496 return result * 0.01; // kA/cm² µm² --> mA
497}
498
499 double ElectricalFem3DSolver::getTotalCurrent(size_t nact) {
500 if (nact >= active.size()) throw BadInput(this->getId(), "wrong active region number");
501 const auto& act = active[nact];
502 // Find the average of the active region
503 size_t level = (act.bottom + act.top) / 2;
504 return integrateCurrent(level, true);
505}
506
507 const LazyData<double> ElectricalFem3DSolver::getVoltage(shared_ptr<const MeshD<3>> dest_mesh, InterpolationMethod method) const {
508 if (!potential) throw NoValue("voltage");
509 this->writelog(LOG_DEBUG, "Getting potential");
510 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
511 if (maskedMesh->full())
512 return interpolate(mesh, potential, dest_mesh, method, geometry);
513 else
514 return interpolate(maskedMesh, potential, dest_mesh, method, geometry);
515}
516
517 const LazyData<Vec<3>> ElectricalFem3DSolver::getCurrentDensity(shared_ptr<const MeshD<3>> dest_mesh, InterpolationMethod method) {
518 if (!potential) throw NoValue("current density");
519 this->writelog(LOG_DEBUG, "Getting current density");
520 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
521 InterpolationFlags flags(geometry, InterpolationFlags::Symmetry::NPP, InterpolationFlags::Symmetry::PNP,
522 InterpolationFlags::Symmetry::PPN);
523 if (maskedMesh->full()) {
524 auto result = interpolate(mesh->getElementMesh(), current, dest_mesh, method, flags);
525 return LazyData<Vec<3>>(result.size(), [this, dest_mesh, result, flags](size_t i) {
526 return this->geometry->getChildBoundingBox().contains(flags.wrap(dest_mesh->at(i))) ? result[i] : Vec<3>(0., 0., 0.);
527 });
528 } else {
529 auto result = interpolate(maskedMesh->getElementMesh(), current, dest_mesh, method, flags);
530 return LazyData<Vec<3>>(result.size(), [result](size_t i) {
531 // Masked mesh always returns NaN outside of itself
532 auto val = result[i];
533 return isnan(val) ? Vec<3>(0., 0., 0.) : val;
534 });
535 }
536}
537
538 const LazyData<double> ElectricalFem3DSolver::getHeatDensity(shared_ptr<const MeshD<3>> dest_mesh, InterpolationMethod method) {
539 if (!potential) throw NoValue("heat density");
540 this->writelog(LOG_DEBUG, "Getting heat density");
541 if (!heat) saveHeatDensity(); // we will compute heats only if they are needed
542 if (method == INTERPOLATION_DEFAULT) method = INTERPOLATION_LINEAR;
543 InterpolationFlags flags(geometry);
544 if (maskedMesh->full()) {
545 auto result = interpolate(mesh->getElementMesh(), heat, dest_mesh, method, flags);
546 return LazyData<double>(result.size(), [this, dest_mesh, result, flags](size_t i) {
547 return this->geometry->getChildBoundingBox().contains(flags.wrap(dest_mesh->at(i))) ? result[i] : 0.;
548 });
549 } else {
550 auto result = interpolate(maskedMesh->getElementMesh(), heat, dest_mesh, method, flags);
551 return LazyData<double>(result.size(), [result](size_t i) {
552 // Masked mesh always returns NaN outside of itself
553 auto val = result[i];
554 return isnan(val) ? 0. : val;
555 });
556 }
557}
558
559 const LazyData<Tensor2<double>> ElectricalFem3DSolver::getConductivity(shared_ptr<const MeshD<3>> dest_mesh,
560 InterpolationMethod /*method*/) {
561 initCalculation();
562 writelog(LOG_DEBUG, "Getting conductivities");
563 loadConductivity();
564 InterpolationFlags flags(geometry);
565 return interpolate(maskedMesh->getElementMesh(), conds, dest_mesh, INTERPOLATION_NEAREST, flags);
566}
567
568 double ElectricalFem3DSolver::getTotalEnergy() {
569 double W = 0.;
570 auto T = inTemperature(maskedMesh->getElementMesh());
571 for (auto el : maskedMesh->elements()) {
572 size_t lll = el.getLoLoLoIndex();
573 size_t llu = el.getLoLoUpIndex();
574 size_t lul = el.getLoUpLoIndex();
575 size_t luu = el.getLoUpUpIndex();
576 size_t ull = el.getUpLoLoIndex();
577 size_t ulu = el.getUpLoUpIndex();
578 size_t uul = el.getUpUpLoIndex();
579 size_t uuu = el.getUpUpUpIndex();
580 double dvx = -0.25e6 *
581 (-potential[lll] - potential[llu] - potential[lul] - potential[luu] + potential[ull] + potential[ulu] +
582 potential[uul] + potential[uuu]) /
583 (el.getUpper0() - el.getLower0()); // 1e6 - from µm to m
584 double dvy = -0.25e6 *
585 (-potential[lll] - potential[llu] + potential[lul] + potential[luu] - potential[ull] - potential[ulu] +
586 potential[uul] + potential[uuu]) /
587 (el.getUpper1() - el.getLower1()); // 1e6 - from µm to m
588 double dvz = -0.25e6 *
589 (-potential[lll] + potential[llu] - potential[lul] + potential[luu] - potential[ull] + potential[ulu] -
590 potential[uul] + potential[uuu]) /
591 (el.getUpper2() - el.getLower2()); // 1e6 - from µm to m
592 double w = this->geometry->getMaterial(el.getMidpoint())->eps(T[el.getIndex()]) * (dvx * dvx + dvy * dvy + dvz * dvz);
593 double d0 = el.getUpper0() - el.getLower0();
594 double d1 = el.getUpper1() - el.getLower1();
595 double d2 = el.getUpper2() - el.getLower2();
596 // TODO add outsides of computational area
597 W += 0.5e-18 * phys::epsilon0 * d0 * d1 * d2 * w; // 1e-18 µm³ -> m³
598 }
599 return W;
600}
601
602 double ElectricalFem3DSolver::getCapacitance() {
603 if (this->voltage_boundary.size() != 2) {
604 throw BadInput(this->getId(), "cannot estimate applied voltage (exactly 2 voltage boundary conditions required)");
605 }
606
607 double U = voltage_boundary[0].value - voltage_boundary[1].value;
608
609 return 2e12 * getTotalEnergy() / (U * U); // 1e12 F -> pF
610}
611
612 double ElectricalFem3DSolver::getTotalHeat() {
613 double W = 0.;
614 if (!heat) saveHeatDensity(); // we will compute heats only if they are needed
615 for (auto el : this->maskedMesh->elements()) {
616 double d0 = el.getUpper0() - el.getLower0();
617 double d1 = el.getUpper1() - el.getLower1();
618 double d2 = el.getUpper2() - el.getLower2();
619 W += 1e-15 * d0 * d1 * d2 * heat[el.getIndex()]; // 1e-15 µm³ -> m³, W -> mW
620 }
621 if (geometry->isSymmetric(Geometry::DIRECTION_LONG)) W *= 2.;
622 if (geometry->isSymmetric(Geometry::DIRECTION_TRAN)) W *= 2.;
623 return W;
624}
625
626}}} // namespace plask::electrical::shockley