PLaSK library
Loading...
Searching...
No Matches
modal-python.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#define PY_ARRAY_UNIQUE_SYMBOL PLASK_OPTICAL_SLAB_ARRAY_API
15
16#include "modal-python.hpp"
17using namespace plask;
18using namespace plask::python;
19
20#include "fourier2d-python.hpp"
21#include "fourier3d-python.hpp"
22#include "besselcyl-python.hpp"
23using namespace plask::optical::modal;
24using namespace plask::optical::modal::python;
25
26#ifndef NDEBUG
27template <typename T>
28struct Matrix_Python {
29 Matrix<T> data;
30 Matrix_Python(const Matrix<T>& data): data(data) {}
31 Matrix_Python(const Matrix_Python<T>& src): data(src.data) {}
32
33 static PyObject* convert(const Matrix<T>& self) {
34 npy_intp dims[2] = { static_cast<npy_intp>(self.rows()), static_cast<npy_intp>(self.cols()) };
35 npy_intp strides[2] = { sizeof(dcomplex), npy_intp(self.rows() * sizeof(dcomplex)) };
36
37 PyObject* arr = PyArray_New(&PyArray_Type, 2, dims, plask::python::detail::typenum<T>(), strides,
38 (void*)self.data(), 0, 0, NULL);
39 if (arr == nullptr) throw plask::CriticalException(u8"cannot create array from matrix");
40 // Make sure the data vector stays alive as long as the array
41 py::object oself {Matrix_Python<T>(self)};
42 py::incref(oself.ptr());
43 PyArray_SetBaseObject((PyArrayObject*)arr, oself.ptr());
44 return arr;
45 }
46};
47
48template <typename T>
49struct Diagonal_Python {
50 MatrixDiagonal<T> data;
51 Diagonal_Python(const MatrixDiagonal<T>& data): data(data) {}
52 Diagonal_Python(const Diagonal_Python<T>& src): data(src.data) {}
53
54 static PyObject* convert(const MatrixDiagonal<T>& self) {
55 npy_intp dims[1] = { static_cast<npy_intp>(self.size()) };
56 npy_intp strides[1] = { sizeof(dcomplex) };
57
58 PyObject* arr = PyArray_New(&PyArray_Type, 1, dims, plask::python::detail::typenum<T>(), strides,
59 (void*)self.data(), 0, 0, NULL);
60 if (arr == nullptr) throw plask::CriticalException(u8"cannot create array from matrix");
61 // Make sure the data vector stays alive as long as the array
62 py::object oself {Diagonal_Python<T>(self)};
63 py::incref(oself.ptr());
64 PyArray_SetBaseObject((PyArrayObject*)arr, oself.ptr());
65 return arr;
66 }
67};
68#endif
69
70
71void* CoeffsArray::convertible(PyObject* obj) {
72 if (PyArray_Check(obj))
73 return obj;
74 return nullptr;
75
76}
77
78void CoeffsArray::construct(PyObject* obj, boost::python::converter::rvalue_from_python_stage1_data* data) {
79 PyArrayObject* arr = (PyArrayObject*)obj;
80 if (PyArray_NDIM(arr) != 1)
81 throw TypeError("only rank 1 arrays allowed");
82 size_t size = PyArray_DIMS(arr)[0];
83
84 if (PyArray_TYPE(arr) != NPY_CDOUBLE || PyArray_STRIDES(arr)[0] != sizeof(dcomplex)) {
85 writelog(LOG_DEBUG, u8"Copying numpy array to make it contiguous");
86 npy_intp sizes[] = { (npy_int)size };
87 npy_intp strides[] = { sizeof(dcomplex) };
88 PyObject* newarr = PyArray_New(&PyArray_Type, 1, sizes,
89 PyArray_TYPE(arr), strides,
90 nullptr, 0, 0, nullptr);
91 PyArray_CopyInto((PyArrayObject*)newarr, arr);
92 arr = (PyArrayObject*)newarr;
93 }
94
95 void* storage = ((boost::python::converter::rvalue_from_python_storage<CoeffsArray>*)data)->storage.bytes;
96 new(storage) CoeffsArray(arr);
97 data->convertible = storage;
98}
99
100
101BOOST_PYTHON_MODULE(modal)
102{
103 plask_import_array();
104
105#ifndef NDEBUG
106 py::class_<Matrix_Python<dcomplex>>("_cmatrix", py::no_init);
107 py::delattr(py::scope(), "_cmatrix");
108 py::to_python_converter<cmatrix, Matrix_Python<dcomplex>>();
109 py::class_<Diagonal_Python<dcomplex>>("_cdiagonal", py::no_init);
110 py::delattr(py::scope(), "_cdiagonal");
111 py::to_python_converter<cdiagonal, Diagonal_Python<dcomplex>>();
112 py::class_<Matrix_Python<double>>("_dmatrix", py::no_init);
113 py::delattr(py::scope(), "_dmatrix");
114 py::to_python_converter<dmatrix, Matrix_Python<double>>();
115 // py::class_<Diagonal_Python<double>>("_ddiagonal", py::no_init);
116 // py::delattr(py::scope(), "_ddiagonal");
117 // py::to_python_converter<MatrixDiagonal<double>, Diagonal_Python<double>>();
118#endif
119
120 py::to_python_converter<Expansion::Component, PythonComponentConventer>();
121 py::converter::registry::push_back(&PythonComponentConventer::convertible, &PythonComponentConventer::construct,
122 py::type_id<Expansion::Component>());
123
124 // py::converter::registry::push_back(&PythonFourierSolver3DWhatConverter::convertible, &PythonFourierSolver3DWhatConverter::construct,
125 // py::type_id<FourierSolver3D::What>());
126
127 py::class_<PmlWrapper, shared_ptr<PmlWrapper>>("PML", "Perfectly matched layer details.", py::no_init)
128 .def("__init__", py::make_constructor(&PmlWrapper::__init__, py::default_call_policies(),
129 (py::arg("factor"), "size", "dist", py::arg("shape")=2)))
130 .add_property("factor", &PmlWrapper::get_factor, &PmlWrapper::set_factor, "PML scaling factor.")
131 .add_property("size", &PmlWrapper::get_size, &PmlWrapper::set_size, u8"PML size.")
132 .add_property("dist", &PmlWrapper::get_dist, &PmlWrapper::set_dist, u8"PML distance from the structure.")
133 .add_property("shape", &PmlWrapper::get_order, &PmlWrapper::set_order, u8"PML shape order (0 → flat, 1 → linearly increasing, 2 → quadratic, etc.).")
134 .def("__str__", &PmlWrapper::__str__)
135 .def("__repr__", &PmlWrapper::__repr__)
136 ;
137
138 py_enum<Transfer::Method>()
139 .value("AUTO", Transfer::METHOD_AUTO)
140 .value("REFLECTION", Transfer::METHOD_REFLECTION_ADMITTANCE)
141 .value("REFLECTION_ADMITTANCE", Transfer::METHOD_REFLECTION_ADMITTANCE)
142 .value("REFLECTION_IMPEDANCE", Transfer::METHOD_REFLECTION_IMPEDANCE)
143 .value("REFLECTION", Transfer::METHOD_REFLECTION_ADMITTANCE)
144 .value("ADMITTANCE", Transfer::METHOD_ADMITTANCE)
145 .value("IMPEDANCE", Transfer::METHOD_IMPEDANCE)
146 ;
147
148 py_enum<Transfer::Determinant>()
149 .value("EIGEN", Transfer::DETERMINANT_EIGENVALUE)
150 .value("EIGENVALUE", Transfer::DETERMINANT_EIGENVALUE)
151 .value("FULL", Transfer::DETERMINANT_FULL)
152 ;
153
154 py_enum<RootDigger::Method>()
155 .value("MULLER", RootDigger::ROOT_MULLER)
156 .value("BROYDEN", RootDigger::ROOT_BROYDEN)
157 .value("BRENT", RootDigger::ROOT_BRENT)
158 ;
159
160 py_enum<typename Transfer::IncidentDirection>()
161 .value("TOP", Transfer::INCIDENCE_TOP)
162 .value("BOTTOM", Transfer::INCIDENCE_BOTTOM)
163 ;
164
165 py_enum<ModalBase::Emission>()
166 .value("UNDEFINED", ModalBase::EMISSION_UNSPECIFIED)
167 .value("TOP", ModalBase::EMISSION_TOP)
168 .value("BOTTOM", ModalBase::EMISSION_BOTTOM)
169 .value("FRONT", ModalBase::EMISSION_FRONT)
170 .value("BACK", ModalBase::EMISSION_BACK)
171 ;
172
173 py::class_<RootDigger::Params, boost::noncopyable>("RootParams", u8"Configuration of the root finding algorithm.", py::no_init)
174 .def_readwrite("method", &RootDigger::Params::method, u8"Root finding method ('muller', 'broyden', or 'brent')")
175 .def_readwrite("tolx", &RootDigger::Params::tolx, u8"Absolute tolerance on the argument.")
176 .def_readwrite("tolf_min", &RootDigger::Params::tolf_min, u8"Sufficient tolerance on the function value.")
177 .def_readwrite("tolf_max", &RootDigger::Params::tolf_max, u8"Required tolerance on the function value.")
178 .def_readwrite("maxiter", &RootDigger::Params::maxiter, u8"Maximum number of iterations.")
179 .def_readwrite("maxstep", &RootDigger::Params::maxstep, u8"Maximum step in one iteration (Broyden method only).")
180 .def_readwrite("alpha", &RootDigger::Params::maxstep, u8"Parameter ensuring sufficient decrease of determinant in each step\n(Broyden method only).")
181 .def_readwrite("lambd", &RootDigger::Params::maxstep, u8"Minimum decrease ratio of one step (Broyden method only).")
182 .def_readwrite("initial_range", &RootDigger::Params::initial_dist, u8"Initial range size (Muller and Brent methods only).")
183 ;
184
185 export_FourierSolver2D();
186 export_FourierSolver3D();
187 export_BesselSolverCyl();
188
189 py::converter::registry::push_back(&CoeffsArray::convertible, &CoeffsArray::construct, boost::python::type_id<CoeffsArray>());
190}