Skip to content

datahandler2

Two-Dimensional Data Handler Implementation for PDE Solvers.

This module implements data handling functionality for 2D PDE problems, focusing on efficient tensor conversions and management of finite element data structures. It provides methods for converting numpy arrays to tensorflow tensors and handling various aspects of the PDE solution process.

The implementation supports
  • Shape function and gradient tensor management
  • Dirichlet boundary data processing
  • Test point generation and handling
  • Sensor data management for inverse problems
  • Bilinear parameter tensor conversion
  • Forcing function data handling
Key classes
  • DataHandler2D: Implementation for 2D PDE data handling
Dependencies
  • tensorflow: For tensor operations
  • numpy: For numerical arrays
  • FESpace2D: For finite element space handling
  • Domain2D: For domain management
Note

The implementation follows FastVPINNs methodology [1] for efficient handling of finite element data structures.

References

[1] FastVPINNs: Tensor-Driven Acceleration of VPINNs for Complex Geometries DOI: https://arxiv.org/abs/2404.12063

DataHandler2D

Bases: DataHandler

Handles data conversion and management for 2D PDE problems.

This class implements the DataHandler interface for 2D problems, providing methods for converting finite element data to tensorflow tensors and managing various aspects of the PDE solution process.

Attributes:

Name Type Description
fespace

Finite element space object for mesh and element info
domain

Domain object for geometric information
dtype

TensorFlow data type for tensor conversion
shape_val_mat_list

Tensor of shape function values Shape: List of matrices of shape (n_test_functions, n_quad_points) with length n_elements
grad_x_mat_list

Tensor of x-derivatives Shape: List of matrices of shape (n_test_functions, n_quad_points) with length n_elements
grad_y_mat_list

Tensor of y-derivatives Shape: List of matrices of shape (n_test_functions, n_quad_points) with length n_elements
x_pde_list

Tensor of quadrature point coordinates
forcing_function_list

Tensor of forcing function values
test_points

Tensor of test point coordinates
Example

fespace = FESpace2D(mesh, elements) domain = Domain2D(bounds) handler = DataHandler2D(fespace, domain, tf.float32) dirichlet_input, dirichlet_vals = handler.get_dirichlet_input() test_points = handler.get_test_points()

Note

All input numpy arrays are assumed to be float64. The class handles conversion to the specified tensorflow dtype (typically float32) for computational efficiency.

Source code in scirex/core/sciml/fastvpinns/data/datahandler2d.py 
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
class DataHandler2D(DataHandler):
    """Handles data conversion and management for 2D PDE problems.

    This class implements the DataHandler interface for 2D problems,
    providing methods for converting finite element data to tensorflow
    tensors and managing various aspects of the PDE solution process.

    Attributes:
        fespace: Finite element space object for mesh and element info
        domain: Domain object for geometric information
        dtype: TensorFlow data type for tensor conversion
        shape_val_mat_list: Tensor of shape function values
            Shape: List of matrices of shape (n_test_functions, n_quad_points) with length n_elements
        grad_x_mat_list: Tensor of x-derivatives
            Shape: List of matrices of shape (n_test_functions, n_quad_points) with length n_elements
        grad_y_mat_list: Tensor of y-derivatives
            Shape: List of matrices of shape (n_test_functions, n_quad_points) with length n_elements
        x_pde_list: Tensor of quadrature point coordinates
        forcing_function_list: Tensor of forcing function values
        test_points: Tensor of test point coordinates

    Example:
        >>> fespace = FESpace2D(mesh, elements)
        >>> domain = Domain2D(bounds)
        >>> handler = DataHandler2D(fespace, domain, tf.float32)
        >>> dirichlet_input, dirichlet_vals = handler.get_dirichlet_input()
        >>> test_points = handler.get_test_points()

    Note:
        All input numpy arrays are assumed to be float64. The class handles
        conversion to the specified tensorflow dtype (typically float32)
        for computational efficiency.
    """

    def __init__(self, fespace, domain, dtype):
        """
        Constructor for the DataHandler2D class

        Args:
            fespace (FESpace2D): The FESpace2D object.
            domain (Domain2D): The Domain2D object.
            dtype (tf.DType): The tensorflow dtype to be used for all the tensors.

        Returns:
            None
        """
        # call the parent class constructor
        super().__init__(fespace=fespace, domain=domain, dtype=dtype)

        self.shape_val_mat_list = []
        self.grad_x_mat_list = []
        self.grad_y_mat_list = []
        self.x_pde_list = []
        self.forcing_function_list = []

        # check if the given dtype is a valid tensorflow dtype
        if not isinstance(self.dtype, tf.DType):
            raise TypeError("The given dtype is not a valid tensorflow dtype")

        for cell_index in range(self.fespace.n_cells):
            shape_val_mat = tf.constant(
                self.fespace.get_shape_function_val(cell_index), dtype=self.dtype
            )
            grad_x_mat = tf.constant(
                self.fespace.get_shape_function_grad_x(cell_index), dtype=self.dtype
            )
            grad_y_mat = tf.constant(
                self.fespace.get_shape_function_grad_y(cell_index), dtype=self.dtype
            )
            x_pde = tf.constant(
                self.fespace.get_quadrature_actual_coordinates(cell_index),
                dtype=self.dtype,
            )
            forcing_function = tf.constant(
                self.fespace.get_forcing_function_values(cell_index), dtype=self.dtype
            )
            self.shape_val_mat_list.append(shape_val_mat)
            self.grad_x_mat_list.append(grad_x_mat)
            self.grad_y_mat_list.append(grad_y_mat)
            self.x_pde_list.append(x_pde)
            self.forcing_function_list.append(forcing_function)

        # now convert all the shapes into 3D tensors for easy multiplication
        # input tensor - x_pde_list
        self.x_pde_list = tf.reshape(self.x_pde_list, [-1, 2])

        self.forcing_function_list = tf.concat(self.forcing_function_list, axis=1)

        self.shape_val_mat_list = tf.stack(self.shape_val_mat_list, axis=0)
        self.grad_x_mat_list = tf.stack(self.grad_x_mat_list, axis=0)
        self.grad_y_mat_list = tf.stack(self.grad_y_mat_list, axis=0)

        # test points
        self.test_points = None

    def get_dirichlet_input(self) -> tuple:
        """
        This function will return the input for the Dirichlet boundary data

        Args:
            None

        Returns:
            The Dirichlet boundary data as a tuple of tensors
        """
        input_dirichlet, actual_dirichlet = (
            self.fespace.generate_dirichlet_boundary_data()
        )

        # convert to tensors
        input_dirichlet = tf.constant(input_dirichlet, dtype=self.dtype)
        actual_dirichlet = tf.constant(actual_dirichlet, dtype=self.dtype)
        actual_dirichlet = tf.reshape(actual_dirichlet, [-1, 1])

        return input_dirichlet, actual_dirichlet

    def get_test_points(self) -> tf.Tensor:
        """
        Get the test points for the given domain.

        Args:
            None

        Returns:
            The test points as a tensor
        """
        self.test_points = self.domain.get_test_points()
        self.test_points = tf.constant(self.test_points, dtype=self.dtype)
        return self.test_points

    def get_bilinear_params_dict_as_tensors(self, function):
        """
        Accepts a function from example file and converts all the values into tensors of the given dtype

        Args:
            function (function): The function from the example file which returns the bilinear parameters dictionary

        Returns:
            The bilinear parameters dictionary with all the values converted to tensors
        """
        # get the dictionary of bilinear parameters
        bilinear_params_dict = function()

        # loop over all keys and convert the values to tensors
        for key in bilinear_params_dict.keys():
            bilinear_params_dict[key] = tf.constant(
                bilinear_params_dict[key], dtype=self.dtype
            )

        return bilinear_params_dict

    # to be used only in inverse problems
    def get_sensor_data(
        self, exact_sol, num_sensor_points, mesh_type, file_name=None
    ) -> tuple:
        """
        Accepts a function from example file and converts all the values into tensors of the given dtype

        Args:
            exact_sol (function): The exact solution function
            num_sensor_points (int): The number of sensor points
            mesh_type (str): The type of mesh
            file_name (str): The file name to save the sensor data

        Returns:
            The sensor data as a tensor

        Raises:
            ValueError: If the mesh type is not internal or external
        """
        print(f"mesh_type = {mesh_type}")
        if mesh_type == "internal":
            # Call the method to get the sensor data
            points, sensor_values = self.fespace.get_sensor_data(
                exact_sol, num_sensor_points
            )
        elif mesh_type == "external":
            # Call the method to get the sensor data
            points, sensor_values = self.fespace.get_sensor_data_external(
                exact_sol, num_sensor_points, file_name
            )
        # convert the points and sensor values into tensors
        points = tf.constant(points, dtype=self.dtype)
        sensor_values = tf.constant(sensor_values, dtype=self.dtype)

        sensor_values = tf.reshape(sensor_values, [-1, 1])
        points = tf.reshape(points, [-1, 2])

        return points, sensor_values

    # get inverse param dict as tensors
    def get_inverse_params(self, inverse_params_dict_function) -> dict:
        """
        Accepts a function from example file and converts all the values into tensors of the given dtype

        Args:
            inverse_params_dict_function (function): The function from the example file which returns the inverse parameters dictionary

        Returns:
            The inverse parameters dictionary with all the values converted to tensors
        """
        # loop over all keys and convert the values to tensors

        inverse_params_dict = inverse_params_dict_function()

        for key in inverse_params_dict.keys():
            inverse_params_dict[key] = tf.constant(
                inverse_params_dict[key], dtype=self.dtype
            )

        return inverse_params_dict

__init__(fespace, domain, dtype)

Constructor for the DataHandler2D class

Parameters:

Name Type Description Default
fespace FESpace2D

The FESpace2D object.

 required
domain Domain2D

The Domain2D object.

 required
dtype DType

The tensorflow dtype to be used for all the tensors.

 required

Returns:

Type Description

None
Source code in scirex/core/sciml/fastvpinns/data/datahandler2d.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def __init__(self, fespace, domain, dtype):
    """
    Constructor for the DataHandler2D class

    Args:
        fespace (FESpace2D): The FESpace2D object.
        domain (Domain2D): The Domain2D object.
        dtype (tf.DType): The tensorflow dtype to be used for all the tensors.

    Returns:
        None
    """
    # call the parent class constructor
    super().__init__(fespace=fespace, domain=domain, dtype=dtype)

    self.shape_val_mat_list = []
    self.grad_x_mat_list = []
    self.grad_y_mat_list = []
    self.x_pde_list = []
    self.forcing_function_list = []

    # check if the given dtype is a valid tensorflow dtype
    if not isinstance(self.dtype, tf.DType):
        raise TypeError("The given dtype is not a valid tensorflow dtype")

    for cell_index in range(self.fespace.n_cells):
        shape_val_mat = tf.constant(
            self.fespace.get_shape_function_val(cell_index), dtype=self.dtype
        )
        grad_x_mat = tf.constant(
            self.fespace.get_shape_function_grad_x(cell_index), dtype=self.dtype
        )
        grad_y_mat = tf.constant(
            self.fespace.get_shape_function_grad_y(cell_index), dtype=self.dtype
        )
        x_pde = tf.constant(
            self.fespace.get_quadrature_actual_coordinates(cell_index),
            dtype=self.dtype,
        )
        forcing_function = tf.constant(
            self.fespace.get_forcing_function_values(cell_index), dtype=self.dtype
        )
        self.shape_val_mat_list.append(shape_val_mat)
        self.grad_x_mat_list.append(grad_x_mat)
        self.grad_y_mat_list.append(grad_y_mat)
        self.x_pde_list.append(x_pde)
        self.forcing_function_list.append(forcing_function)

    # now convert all the shapes into 3D tensors for easy multiplication
    # input tensor - x_pde_list
    self.x_pde_list = tf.reshape(self.x_pde_list, [-1, 2])

    self.forcing_function_list = tf.concat(self.forcing_function_list, axis=1)

    self.shape_val_mat_list = tf.stack(self.shape_val_mat_list, axis=0)
    self.grad_x_mat_list = tf.stack(self.grad_x_mat_list, axis=0)
    self.grad_y_mat_list = tf.stack(self.grad_y_mat_list, axis=0)

    # test points
    self.test_points = None

get_bilinear_params_dict_as_tensors(function)

Accepts a function from example file and converts all the values into tensors of the given dtype

Parameters:

Name Type Description Default
function function

The function from the example file which returns the bilinear parameters dictionary

 required

Returns:

Type Description

The bilinear parameters dictionary with all the values converted to tensors
Source code in scirex/core/sciml/fastvpinns/data/datahandler2d.py 
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
def get_bilinear_params_dict_as_tensors(self, function):
    """
    Accepts a function from example file and converts all the values into tensors of the given dtype

    Args:
        function (function): The function from the example file which returns the bilinear parameters dictionary

    Returns:
        The bilinear parameters dictionary with all the values converted to tensors
    """
    # get the dictionary of bilinear parameters
    bilinear_params_dict = function()

    # loop over all keys and convert the values to tensors
    for key in bilinear_params_dict.keys():
        bilinear_params_dict[key] = tf.constant(
            bilinear_params_dict[key], dtype=self.dtype
        )

    return bilinear_params_dict

get_dirichlet_input()

This function will return the input for the Dirichlet boundary data

Returns:

Type Description
tuple

The Dirichlet boundary data as a tuple of tensors
Source code in scirex/core/sciml/fastvpinns/data/datahandler2d.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def get_dirichlet_input(self) -> tuple:
    """
    This function will return the input for the Dirichlet boundary data

    Args:
        None

    Returns:
        The Dirichlet boundary data as a tuple of tensors
    """
    input_dirichlet, actual_dirichlet = (
        self.fespace.generate_dirichlet_boundary_data()
    )

    # convert to tensors
    input_dirichlet = tf.constant(input_dirichlet, dtype=self.dtype)
    actual_dirichlet = tf.constant(actual_dirichlet, dtype=self.dtype)
    actual_dirichlet = tf.reshape(actual_dirichlet, [-1, 1])

    return input_dirichlet, actual_dirichlet

get_inverse_params(inverse_params_dict_function)

Accepts a function from example file and converts all the values into tensors of the given dtype

Parameters:

Name Type Description Default
inverse_params_dict_function function

The function from the example file which returns the inverse parameters dictionary

 required

Returns:

Type Description
dict

The inverse parameters dictionary with all the values converted to tensors
Source code in scirex/core/sciml/fastvpinns/data/datahandler2d.py 
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
def get_inverse_params(self, inverse_params_dict_function) -> dict:
    """
    Accepts a function from example file and converts all the values into tensors of the given dtype

    Args:
        inverse_params_dict_function (function): The function from the example file which returns the inverse parameters dictionary

    Returns:
        The inverse parameters dictionary with all the values converted to tensors
    """
    # loop over all keys and convert the values to tensors

    inverse_params_dict = inverse_params_dict_function()

    for key in inverse_params_dict.keys():
        inverse_params_dict[key] = tf.constant(
            inverse_params_dict[key], dtype=self.dtype
        )

    return inverse_params_dict

get_sensor_data(exact_sol, num_sensor_points, mesh_type, file_name=None)

Accepts a function from example file and converts all the values into tensors of the given dtype

Parameters:

Name Type Description Default
exact_sol function

The exact solution function

 required
num_sensor_points int

The number of sensor points

 required
mesh_type str

The type of mesh

 required
file_name str

The file name to save the sensor data

 None

Returns:

Type Description
tuple

The sensor data as a tensor

Raises:

Type Description
ValueError

If the mesh type is not internal or external
Source code in scirex/core/sciml/fastvpinns/data/datahandler2d.py 
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
def get_sensor_data(
    self, exact_sol, num_sensor_points, mesh_type, file_name=None
) -> tuple:
    """
    Accepts a function from example file and converts all the values into tensors of the given dtype

    Args:
        exact_sol (function): The exact solution function
        num_sensor_points (int): The number of sensor points
        mesh_type (str): The type of mesh
        file_name (str): The file name to save the sensor data

    Returns:
        The sensor data as a tensor

    Raises:
        ValueError: If the mesh type is not internal or external
    """
    print(f"mesh_type = {mesh_type}")
    if mesh_type == "internal":
        # Call the method to get the sensor data
        points, sensor_values = self.fespace.get_sensor_data(
            exact_sol, num_sensor_points
        )
    elif mesh_type == "external":
        # Call the method to get the sensor data
        points, sensor_values = self.fespace.get_sensor_data_external(
            exact_sol, num_sensor_points, file_name
        )
    # convert the points and sensor values into tensors
    points = tf.constant(points, dtype=self.dtype)
    sensor_values = tf.constant(sensor_values, dtype=self.dtype)

    sensor_values = tf.reshape(sensor_values, [-1, 1])
    points = tf.reshape(points, [-1, 2])

    return points, sensor_values

get_test_points()

Get the test points for the given domain.

Returns:

Type Description
Tensor

The test points as a tensor
Source code in scirex/core/sciml/fastvpinns/data/datahandler2d.py 
181
182
183
184
185
186
187
188
189
190
191
192
193
def get_test_points(self) -> tf.Tensor:
    """
    Get the test points for the given domain.

    Args:
        None

    Returns:
        The test points as a tensor
    """
    self.test_points = self.domain.get_test_points()
    self.test_points = tf.constant(self.test_points, dtype=self.dtype)
    return self.test_points

FE2D_Cell

A class for managing finite element computations at the cell level.

This class handles the storage and computation of finite element values, including basis functions, quadrature rules, and transformations for a single cell in a 2D mesh.

Attributes:

Name Type Description
cell_coordinates ndarray

Physical coordinates of the cell vertices
cell_type str

Type of the cell (e.g., 'quad', 'triangle')
fe_order int

Order of the finite element approximation
fe_type str

Type of finite element basis
quad_order int

Order of quadrature rule
quad_type str

Type of quadrature formula
fe_transformation str

Type of geometric transformation
forcing_function callable

Source term function
basis_function BasisFunction2D

Basis function implementation
quad_xi ndarray

Xi coordinates of quadrature points
quad_eta ndarray

Eta coordinates of quadrature points
quad_weight ndarray

Quadrature weights
jacobian ndarray

Transformation Jacobian
basis_at_quad ndarray

Basis values at quadrature points
basis_gradx_at_quad ndarray

X-derivatives at quadrature points
basis_grady_at_quad ndarray

Y-derivatives at quadrature points
quad_actual_coordinates ndarray

Physical quadrature point coordinates
Example

coords = np.array([[0,0], [1,0], [1,1], [0,1]]) cell = FE2D_Cell( ... cell_coordinates=coords, ... cell_type='quad', ... fe_order=2, ... fe_type='lagrange', ... quad_order=3, ... quad_type='gauss', ... fe_transformation_type='bilinear', ... forcing_function=lambda x, y: x*y ... ) cell.basis_at_quad # Get basis values at quadrature points

Notes
  • All gradient and derivative values are stored in the reference domain
  • Jacobian and quadrature weights are combined for efficiency
  • Forcing function values are typically computed in the fespace class
  • Supports multiple types of transformations and element types
Source code in scirex/core/sciml/fe/fe2d_cell.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
class FE2D_Cell:
    """
    A class for managing finite element computations at the cell level.

    This class handles the storage and computation of finite element values,
    including basis functions, quadrature rules, and transformations for a
    single cell in a 2D mesh.

    Attributes:
        cell_coordinates (np.ndarray): Physical coordinates of the cell vertices
        cell_type (str): Type of the cell (e.g., 'quad', 'triangle')
        fe_order (int): Order of the finite element approximation
        fe_type (str): Type of finite element basis
        quad_order (int): Order of quadrature rule
        quad_type (str): Type of quadrature formula
        fe_transformation (str): Type of geometric transformation
        forcing_function (callable): Source term function
        basis_function (BasisFunction2D): Basis function implementation
        quad_xi (np.ndarray): Xi coordinates of quadrature points
        quad_eta (np.ndarray): Eta coordinates of quadrature points
        quad_weight (np.ndarray): Quadrature weights
        jacobian (np.ndarray): Transformation Jacobian
        basis_at_quad (np.ndarray): Basis values at quadrature points
        basis_gradx_at_quad (np.ndarray): X-derivatives at quadrature points
        basis_grady_at_quad (np.ndarray): Y-derivatives at quadrature points
        quad_actual_coordinates (np.ndarray): Physical quadrature point coordinates

    Example:
        >>> coords = np.array([[0,0], [1,0], [1,1], [0,1]])
        >>> cell = FE2D_Cell(
        ...     cell_coordinates=coords,
        ...     cell_type='quad',
        ...     fe_order=2,
        ...     fe_type='lagrange',
        ...     quad_order=3,
        ...     quad_type='gauss',
        ...     fe_transformation_type='bilinear',
        ...     forcing_function=lambda x, y: x*y
        ... )
        >>> cell.basis_at_quad  # Get basis values at quadrature points

    Notes:
        - All gradient and derivative values are stored in the reference domain
        - Jacobian and quadrature weights are combined for efficiency
        - Forcing function values are typically computed in the fespace class
        - Supports multiple types of transformations and element types
    """

    def __init__(
        self,
        cell_coordinates: np.ndarray,
        cell_type: str,
        fe_order: int,
        fe_type: str,
        quad_order: int,
        quad_type: str,
        fe_transformation_type: str,
        forcing_function,
    ):
        """
        Constructor for the FE2D_Cell class.

        Args:
            cell_coordinates (np.ndarray): Physical coordinates of the cell vertices
            cell_type (str): Type of the cell (e.g., 'quad', 'triangle')
            fe_order (int): Order of the finite element approximation
            fe_type (str): Type of finite element basis
            quad_order (int): Order of quadrature rule
            quad_type (str): Type of quadrature formula
            fe_transformation_type (str): Type of geometric transformation
            forcing_function (callable): Source term function

        Returns:
            None
        """
        self.cell_coordinates = cell_coordinates
        self.cell_type = cell_type
        self.fe_order = fe_order
        self.fe_type = fe_type
        self.quad_order = quad_order
        self.quad_type = quad_type
        self.fe_transformation = fe_transformation_type
        self.forcing_function = forcing_function

        # Basis function Class
        self.basis_function = None

        # Quadrature Values
        self.quad_xi = None
        self.quad_eta = None
        self.quad_weight = None
        self.jacobian = None
        self.mult = None

        # FE Values
        self.basis_at_quad = None
        self.basis_gradx_at_quad = None
        self.basis_grady_at_quad = None
        self.basis_gradxy_at_quad = None
        self.basis_gradxx_at_quad = None
        self.basis_gradyy_at_quad = None

        # Quadrature Coordinates
        self.quad_actual_coordinates = None

        # Forcing function values at the quadrature points
        self.forcing_at_quad = None

        # FE Transformation Class
        self.fetransformation = None

        # get instance of the FE_setup class
        self.fe_setup = FE2DSetupMain(
            cell_type=self.cell_type,
            fe_order=self.fe_order,
            fe_type=self.fe_type,
            quad_order=self.quad_order,
            quad_type=self.quad_type,
        )

        # Call the function to assign the basis function
        self.assign_basis_function()

        # Assign the quadrature points and weights
        self.assign_quadrature()

        # Assign the FE Transformation
        self.assign_fe_transformation()

        # calculate mult -> quadrature weights * Jacobian
        self.assign_quad_weights_and_jacobian()

        # Calculate the basis function values at the quadrature points
        self.assign_basis_values_at_quadrature_points()

        # calculate the actual coordinates of the quadrature points
        self.assign_quadrature_coordinates()

        # Calculate the forcing function values at the actual quadrature points
        # NOTE : The function is just for printing the shape of the force matrix, the
        # actual calculation is performed on the fespace class
        self.assign_forcing_term(self.forcing_function)

        # # print the values
        # print("============================================================================")
        # print("Cell Co-ord : ", self.cell_coordinates)
        # print("Basis function values at the quadrature points: \n", self.basis_at_quad / self.mult)
        # print("Basis function gradx at the quadrature points: \n", self.basis_gradx_at_quad)
        # print("Basis function grady at the quadrature points: \n", self.basis_grady_at_quad)
        # print("Forcing function values at the quadrature points: \n", self.forcing_at_quad)

        # grad_x = np.array([5,6,7,8])
        # grad_y = np.array([1,2,3,4])

        # pde = np.matmul(self.basis_gradx_at_quad, grad_x.reshape(-1,1)) + np.matmul(self.basis_grady_at_quad, grad_y.reshape(-1,1))
        # print("PDE values at the quadrature points: \n", pde)

    def assign_basis_function(self) -> BasisFunction2D:
        """
        Assigns the basis function class based on the cell type and the FE order.

        Args:
            None

        Returns:
            BasisFunction2D: The basis function class for the given cell type and FE order.
        """
        self.basis_function = self.fe_setup.assign_basis_function()

    def assign_quadrature(self) -> None:
        """
        Assigns the quadrature points and weights based on the cell type and the quadrature order.

        Args:
            None

        Returns:
            None
        """
        self.quad_weight, self.quad_xi, self.quad_eta = (
            self.fe_setup.assign_quadrature_rules()
        )

    def assign_fe_transformation(self) -> None:
        """
        Assigns the FE Transformation class based on the cell type and the FE order.

        This method assigns the appropriate FE Transformation class based on the cell type and the FE order.
        It sets the cell coordinates for the FE Transformation and obtains the Jacobian of the transformation.

        Args:
            None

        Returns:
            None
        """
        self.fetransformation = self.fe_setup.assign_fe_transformation(
            self.fe_transformation, self.cell_coordinates
        )
        # Sets cell co-ordinates for the FE Transformation
        self.fetransformation.set_cell()

        # obtains the Jacobian of the transformation
        self.jacobian = self.fetransformation.get_jacobian(
            self.quad_xi, self.quad_eta
        ).reshape(-1, 1)

    def assign_basis_values_at_quadrature_points(self) -> None:
        """
        Assigns the basis function values at the quadrature points.

        This method calculates the values of the basis functions and their gradients at the quadrature points.
        The basis function values are stored in `self.basis_at_quad`, while the gradients are stored in
        `self.basis_gradx_at_quad`, `self.basis_grady_at_quad`, `self.basis_gradxy_at_quad`,
        `self.basis_gradxx_at_quad`, and `self.basis_gradyy_at_quad`.

        Args:
            None

        Returns:
            None
        """
        self.basis_at_quad = []
        self.basis_gradx_at_quad = []
        self.basis_grady_at_quad = []
        self.basis_gradxy_at_quad = []
        self.basis_gradxx_at_quad = []
        self.basis_gradyy_at_quad = []

        self.basis_at_quad = self.basis_function.value(self.quad_xi, self.quad_eta)

        # For Gradients we need to perform a transformation to the original cell
        grad_x_ref = self.basis_function.gradx(self.quad_xi, self.quad_eta)
        grad_y_ref = self.basis_function.grady(self.quad_xi, self.quad_eta)

        grad_x_orig, grad_y_orig = self.fetransformation.get_orig_from_ref_derivative(
            grad_x_ref, grad_y_ref, self.quad_xi, self.quad_eta
        )

        self.basis_gradx_at_quad = grad_x_orig
        self.basis_grady_at_quad = grad_y_orig

        self.basis_gradx_at_quad_ref = grad_x_ref
        self.basis_grady_at_quad_ref = grad_y_ref

        # get the double derivatives of the basis functions ( ref co-ordinates )
        grad_xx_ref = self.basis_function.gradxx(self.quad_xi, self.quad_eta)
        grad_xy_ref = self.basis_function.gradxy(self.quad_xi, self.quad_eta)
        grad_yy_ref = self.basis_function.gradyy(self.quad_xi, self.quad_eta)

        # get the double derivatives of the basis functions ( orig co-ordinates )
        grad_xx_orig, grad_xy_orig, grad_yy_orig = (
            self.fetransformation.get_orig_from_ref_second_derivative(
                grad_xx_ref, grad_xy_ref, grad_yy_ref, self.quad_xi, self.quad_eta
            )
        )

        # = the value
        self.basis_gradxy_at_quad = grad_xy_orig
        self.basis_gradxx_at_quad = grad_xx_orig
        self.basis_gradyy_at_quad = grad_yy_orig

        # Multiply each row with the quadrature weights
        # Basis at Quad - n_test * N_quad
        self.basis_at_quad = self.basis_at_quad * self.mult
        self.basis_gradx_at_quad = self.basis_gradx_at_quad * self.mult
        self.basis_grady_at_quad = self.basis_grady_at_quad * self.mult
        self.basis_gradxy_at_quad = self.basis_gradxy_at_quad * self.mult
        self.basis_gradxx_at_quad = self.basis_gradxx_at_quad * self.mult
        self.basis_gradyy_at_quad = self.basis_gradyy_at_quad * self.mult

    def assign_quad_weights_and_jacobian(self) -> None:
        """
        Assigns the quadrature weights and the Jacobian of the transformation.

        This method calculates and assigns the quadrature weights and the Jacobian of the transformation
        for the current cell. The quadrature weights are multiplied by the flattened Jacobian array
        and stored in the `mult` attribute of the class.

        Args:
            None

        Returns:
            None
        """
        self.mult = self.quad_weight * self.jacobian.flatten()

    def assign_quadrature_coordinates(self) -> None:
        """
        Assigns the actual coordinates of the quadrature points.

        This method calculates the actual coordinates of the quadrature points based on the given Xi and Eta values.
        The Xi and Eta values are obtained from the `quad_xi` and `quad_eta` attributes of the class.
        The calculated coordinates are stored in the `quad_actual_coordinates` attribute as a NumPy array.

        Args:
            None

        Returns:
            None
        """
        actual_co_ord = []
        for xi, eta in zip(self.quad_xi, self.quad_eta):
            actual_co_ord.append(self.fetransformation.get_original_from_ref(xi, eta))

        self.quad_actual_coordinates = np.array(actual_co_ord)

    def assign_forcing_term(self, forcing_function) -> None:
        """
        Assigns the forcing function values at the quadrature points.

        This function computes the values of the forcing function at the quadrature points
        and assigns them to the `forcing_at_quad` attribute of the FE2D_Cell object.

        Args:
            forcing_function (callable): The forcing function to be integrated

        Returns:
            None

        Notes:
            - The final shape of `forcing_at_quad` will be N_shape_functions x 1.
            - This function is for backward compatibility with old code and currently assigns
              the values as zeros. The actual calculation is performed in the fespace class.
        """
        # get number of shape functions
        n_shape_functions = self.basis_function.num_shape_functions

        # Loop over all the basis functions and compute the integral
        f_integral = np.zeros((n_shape_functions, 1), dtype=np.float64)

        # The above code is for backward compatibility with old code. this function will just assign the values as zeros
        # the actual calculation is performed in the fespace class

        # for i in range(n_shape_functions):
        #     val = 0
        #     for q in range(self.basis_at_quad.shape[1]):
        #         x = self.quad_actual_coordinates[q, 0]
        #         y = self.quad_actual_coordinates[q, 1]
        #         # print("f_values[q] = ",f_values[q])

        #         # the JAcobian and the quadrature weights are pre multiplied to the basis functions
        #         val +=  ( self.basis_at_quad[i, q] ) * self.forcing_function(x, y)
        #         # print("val = ", val)

        #     f_integral[i] = val

        self.forcing_at_quad = f_integral

__init__(cell_coordinates, cell_type, fe_order, fe_type, quad_order, quad_type, fe_transformation_type, forcing_function)

Constructor for the FE2D_Cell class.

Parameters:

Name Type Description Default
cell_coordinates ndarray

Physical coordinates of the cell vertices

 required
cell_type str

Type of the cell (e.g., 'quad', 'triangle')

 required
fe_order int

Order of the finite element approximation

 required
fe_type str

Type of finite element basis

 required
quad_order int

Order of quadrature rule

 required
quad_type str

Type of quadrature formula

 required
fe_transformation_type str

Type of geometric transformation

 required
forcing_function callable

Source term function

 required

Returns:

Type Description

None
Source code in scirex/core/sciml/fe/fe2d_cell.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
def __init__(
    self,
    cell_coordinates: np.ndarray,
    cell_type: str,
    fe_order: int,
    fe_type: str,
    quad_order: int,
    quad_type: str,
    fe_transformation_type: str,
    forcing_function,
):
    """
    Constructor for the FE2D_Cell class.

    Args:
        cell_coordinates (np.ndarray): Physical coordinates of the cell vertices
        cell_type (str): Type of the cell (e.g., 'quad', 'triangle')
        fe_order (int): Order of the finite element approximation
        fe_type (str): Type of finite element basis
        quad_order (int): Order of quadrature rule
        quad_type (str): Type of quadrature formula
        fe_transformation_type (str): Type of geometric transformation
        forcing_function (callable): Source term function

    Returns:
        None
    """
    self.cell_coordinates = cell_coordinates
    self.cell_type = cell_type
    self.fe_order = fe_order
    self.fe_type = fe_type
    self.quad_order = quad_order
    self.quad_type = quad_type
    self.fe_transformation = fe_transformation_type
    self.forcing_function = forcing_function

    # Basis function Class
    self.basis_function = None

    # Quadrature Values
    self.quad_xi = None
    self.quad_eta = None
    self.quad_weight = None
    self.jacobian = None
    self.mult = None

    # FE Values
    self.basis_at_quad = None
    self.basis_gradx_at_quad = None
    self.basis_grady_at_quad = None
    self.basis_gradxy_at_quad = None
    self.basis_gradxx_at_quad = None
    self.basis_gradyy_at_quad = None

    # Quadrature Coordinates
    self.quad_actual_coordinates = None

    # Forcing function values at the quadrature points
    self.forcing_at_quad = None

    # FE Transformation Class
    self.fetransformation = None

    # get instance of the FE_setup class
    self.fe_setup = FE2DSetupMain(
        cell_type=self.cell_type,
        fe_order=self.fe_order,
        fe_type=self.fe_type,
        quad_order=self.quad_order,
        quad_type=self.quad_type,
    )

    # Call the function to assign the basis function
    self.assign_basis_function()

    # Assign the quadrature points and weights
    self.assign_quadrature()

    # Assign the FE Transformation
    self.assign_fe_transformation()

    # calculate mult -> quadrature weights * Jacobian
    self.assign_quad_weights_and_jacobian()

    # Calculate the basis function values at the quadrature points
    self.assign_basis_values_at_quadrature_points()

    # calculate the actual coordinates of the quadrature points
    self.assign_quadrature_coordinates()

    # Calculate the forcing function values at the actual quadrature points
    # NOTE : The function is just for printing the shape of the force matrix, the
    # actual calculation is performed on the fespace class
    self.assign_forcing_term(self.forcing_function)

assign_basis_function()

Assigns the basis function class based on the cell type and the FE order.

Returns:

Name Type Description
BasisFunction2D BasisFunction2D

The basis function class for the given cell type and FE order.
Source code in scirex/core/sciml/fe/fe2d_cell.py 
230
231
232
233
234
235
236
237
238
239
240
def assign_basis_function(self) -> BasisFunction2D:
    """
    Assigns the basis function class based on the cell type and the FE order.

    Args:
        None

    Returns:
        BasisFunction2D: The basis function class for the given cell type and FE order.
    """
    self.basis_function = self.fe_setup.assign_basis_function()

assign_basis_values_at_quadrature_points()

Assigns the basis function values at the quadrature points.

This method calculates the values of the basis functions and their gradients at the quadrature points. The basis function values are stored in self.basis_at_quad, while the gradients are stored in self.basis_gradx_at_quad, self.basis_grady_at_quad, self.basis_gradxy_at_quad, self.basis_gradxx_at_quad, and self.basis_gradyy_at_quad.

Returns:

Type Description
None

None
Source code in scirex/core/sciml/fe/fe2d_cell.py 
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
def assign_basis_values_at_quadrature_points(self) -> None:
    """
    Assigns the basis function values at the quadrature points.

    This method calculates the values of the basis functions and their gradients at the quadrature points.
    The basis function values are stored in `self.basis_at_quad`, while the gradients are stored in
    `self.basis_gradx_at_quad`, `self.basis_grady_at_quad`, `self.basis_gradxy_at_quad`,
    `self.basis_gradxx_at_quad`, and `self.basis_gradyy_at_quad`.

    Args:
        None

    Returns:
        None
    """
    self.basis_at_quad = []
    self.basis_gradx_at_quad = []
    self.basis_grady_at_quad = []
    self.basis_gradxy_at_quad = []
    self.basis_gradxx_at_quad = []
    self.basis_gradyy_at_quad = []

    self.basis_at_quad = self.basis_function.value(self.quad_xi, self.quad_eta)

    # For Gradients we need to perform a transformation to the original cell
    grad_x_ref = self.basis_function.gradx(self.quad_xi, self.quad_eta)
    grad_y_ref = self.basis_function.grady(self.quad_xi, self.quad_eta)

    grad_x_orig, grad_y_orig = self.fetransformation.get_orig_from_ref_derivative(
        grad_x_ref, grad_y_ref, self.quad_xi, self.quad_eta
    )

    self.basis_gradx_at_quad = grad_x_orig
    self.basis_grady_at_quad = grad_y_orig

    self.basis_gradx_at_quad_ref = grad_x_ref
    self.basis_grady_at_quad_ref = grad_y_ref

    # get the double derivatives of the basis functions ( ref co-ordinates )
    grad_xx_ref = self.basis_function.gradxx(self.quad_xi, self.quad_eta)
    grad_xy_ref = self.basis_function.gradxy(self.quad_xi, self.quad_eta)
    grad_yy_ref = self.basis_function.gradyy(self.quad_xi, self.quad_eta)

    # get the double derivatives of the basis functions ( orig co-ordinates )
    grad_xx_orig, grad_xy_orig, grad_yy_orig = (
        self.fetransformation.get_orig_from_ref_second_derivative(
            grad_xx_ref, grad_xy_ref, grad_yy_ref, self.quad_xi, self.quad_eta
        )
    )

    # = the value
    self.basis_gradxy_at_quad = grad_xy_orig
    self.basis_gradxx_at_quad = grad_xx_orig
    self.basis_gradyy_at_quad = grad_yy_orig

    # Multiply each row with the quadrature weights
    # Basis at Quad - n_test * N_quad
    self.basis_at_quad = self.basis_at_quad * self.mult
    self.basis_gradx_at_quad = self.basis_gradx_at_quad * self.mult
    self.basis_grady_at_quad = self.basis_grady_at_quad * self.mult
    self.basis_gradxy_at_quad = self.basis_gradxy_at_quad * self.mult
    self.basis_gradxx_at_quad = self.basis_gradxx_at_quad * self.mult
    self.basis_gradyy_at_quad = self.basis_gradyy_at_quad * self.mult

assign_fe_transformation()

Assigns the FE Transformation class based on the cell type and the FE order.

This method assigns the appropriate FE Transformation class based on the cell type and the FE order. It sets the cell coordinates for the FE Transformation and obtains the Jacobian of the transformation.

Returns:

Type Description
None

None
Source code in scirex/core/sciml/fe/fe2d_cell.py 
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
def assign_fe_transformation(self) -> None:
    """
    Assigns the FE Transformation class based on the cell type and the FE order.

    This method assigns the appropriate FE Transformation class based on the cell type and the FE order.
    It sets the cell coordinates for the FE Transformation and obtains the Jacobian of the transformation.

    Args:
        None

    Returns:
        None
    """
    self.fetransformation = self.fe_setup.assign_fe_transformation(
        self.fe_transformation, self.cell_coordinates
    )
    # Sets cell co-ordinates for the FE Transformation
    self.fetransformation.set_cell()

    # obtains the Jacobian of the transformation
    self.jacobian = self.fetransformation.get_jacobian(
        self.quad_xi, self.quad_eta
    ).reshape(-1, 1)

assign_forcing_term(forcing_function)

Assigns the forcing function values at the quadrature points.

This function computes the values of the forcing function at the quadrature points and assigns them to the forcing_at_quad attribute of the FE2D_Cell object.

Parameters:

Name Type Description Default
forcing_function callable

The forcing function to be integrated

 required

Returns:

Type Description
None

None
Notes
  • The final shape of forcing_at_quad will be N_shape_functions x 1.
  • This function is for backward compatibility with old code and currently assigns the values as zeros. The actual calculation is performed in the fespace class.
Source code in scirex/core/sciml/fe/fe2d_cell.py 
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
def assign_forcing_term(self, forcing_function) -> None:
    """
    Assigns the forcing function values at the quadrature points.

    This function computes the values of the forcing function at the quadrature points
    and assigns them to the `forcing_at_quad` attribute of the FE2D_Cell object.

    Args:
        forcing_function (callable): The forcing function to be integrated

    Returns:
        None

    Notes:
        - The final shape of `forcing_at_quad` will be N_shape_functions x 1.
        - This function is for backward compatibility with old code and currently assigns
          the values as zeros. The actual calculation is performed in the fespace class.
    """
    # get number of shape functions
    n_shape_functions = self.basis_function.num_shape_functions

    # Loop over all the basis functions and compute the integral
    f_integral = np.zeros((n_shape_functions, 1), dtype=np.float64)

    # The above code is for backward compatibility with old code. this function will just assign the values as zeros
    # the actual calculation is performed in the fespace class

    # for i in range(n_shape_functions):
    #     val = 0
    #     for q in range(self.basis_at_quad.shape[1]):
    #         x = self.quad_actual_coordinates[q, 0]
    #         y = self.quad_actual_coordinates[q, 1]
    #         # print("f_values[q] = ",f_values[q])

    #         # the JAcobian and the quadrature weights are pre multiplied to the basis functions
    #         val +=  ( self.basis_at_quad[i, q] ) * self.forcing_function(x, y)
    #         # print("val = ", val)

    #     f_integral[i] = val

    self.forcing_at_quad = f_integral

assign_quad_weights_and_jacobian()

Assigns the quadrature weights and the Jacobian of the transformation.

This method calculates and assigns the quadrature weights and the Jacobian of the transformation for the current cell. The quadrature weights are multiplied by the flattened Jacobian array and stored in the mult attribute of the class.

Returns:

Type Description
None

None
Source code in scirex/core/sciml/fe/fe2d_cell.py 
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
def assign_quad_weights_and_jacobian(self) -> None:
    """
    Assigns the quadrature weights and the Jacobian of the transformation.

    This method calculates and assigns the quadrature weights and the Jacobian of the transformation
    for the current cell. The quadrature weights are multiplied by the flattened Jacobian array
    and stored in the `mult` attribute of the class.

    Args:
        None

    Returns:
        None
    """
    self.mult = self.quad_weight * self.jacobian.flatten()

assign_quadrature()

Assigns the quadrature points and weights based on the cell type and the quadrature order.

Returns:

Type Description
None

None
Source code in scirex/core/sciml/fe/fe2d_cell.py 
242
243
244
245
246
247
248
249
250
251
252
253
254
def assign_quadrature(self) -> None:
    """
    Assigns the quadrature points and weights based on the cell type and the quadrature order.

    Args:
        None

    Returns:
        None
    """
    self.quad_weight, self.quad_xi, self.quad_eta = (
        self.fe_setup.assign_quadrature_rules()
    )

assign_quadrature_coordinates()

Assigns the actual coordinates of the quadrature points.

This method calculates the actual coordinates of the quadrature points based on the given Xi and Eta values. The Xi and Eta values are obtained from the quad_xi and quad_eta attributes of the class. The calculated coordinates are stored in the quad_actual_coordinates attribute as a NumPy array.

Returns:

Type Description
None

None
Source code in scirex/core/sciml/fe/fe2d_cell.py 
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
def assign_quadrature_coordinates(self) -> None:
    """
    Assigns the actual coordinates of the quadrature points.

    This method calculates the actual coordinates of the quadrature points based on the given Xi and Eta values.
    The Xi and Eta values are obtained from the `quad_xi` and `quad_eta` attributes of the class.
    The calculated coordinates are stored in the `quad_actual_coordinates` attribute as a NumPy array.

    Args:
        None

    Returns:
        None
    """
    actual_co_ord = []
    for xi, eta in zip(self.quad_xi, self.quad_eta):
        actual_co_ord.append(self.fetransformation.get_original_from_ref(xi, eta))

    self.quad_actual_coordinates = np.array(actual_co_ord)

Fespace

Abstract base class defining the interface for finite element spaces.

This class provides the foundation for implementing finite element spaces, including mesh handling, element operations, and solution computations.

Attributes:

Name Type Description
mesh

Mesh object containing geometric information
cells ndarray

Array of cell indices
boundary_points dict

Dictionary of boundary point information
cell_type str

Type of finite element cell
fe_order int

Order of finite element approximation
fe_type str

Type of finite element basis
quad_order int

Order of quadrature rule
quad_type str

Type of quadrature formula
fe_transformation_type str

Type of geometric transformation
bound_function_dict dict

Dictionary of boundary condition functions
bound_condition_dict dict

Dictionary of boundary condition types
forcing_function callable

Source term function
output_path str

Path for output files
Example

class MyFespace(Fespace): ... def set_finite_elements(self): ... # Implementation ... pass ... def generate_dirichlet_boundary_data(self): ... # Implementation ... pass ... # Implement other abstract methods

Notes
  • All coordinate transformations must be implemented
  • Shape function values and gradients are available in both reference and physical spaces
  • Supports both internal and external sensor data for inverse problems
  • Boundary conditions must be properly specified through the boundary dictionaries
Source code in scirex/core/sciml/fe/fespace.py 
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
class Fespace:
    """
    Abstract base class defining the interface for finite element spaces.

    This class provides the foundation for implementing finite element spaces,
    including mesh handling, element operations, and solution computations.

    Attributes:
        mesh: Mesh object containing geometric information
        cells (ndarray): Array of cell indices
        boundary_points (dict): Dictionary of boundary point information
        cell_type (str): Type of finite element cell
        fe_order (int): Order of finite element approximation
        fe_type (str): Type of finite element basis
        quad_order (int): Order of quadrature rule
        quad_type (str): Type of quadrature formula
        fe_transformation_type (str): Type of geometric transformation
        bound_function_dict (dict): Dictionary of boundary condition functions
        bound_condition_dict (dict): Dictionary of boundary condition types
        forcing_function (callable): Source term function
        output_path (str): Path for output files

    Example:
        >>> class MyFespace(Fespace):
        ...     def set_finite_elements(self):
        ...         # Implementation
        ...         pass
        ...     def generate_dirichlet_boundary_data(self):
        ...         # Implementation
        ...         pass
        ...     # Implement other abstract methods

    Notes:
        - All coordinate transformations must be implemented
        - Shape function values and gradients are available in both
        reference and physical spaces
        - Supports both internal and external sensor data for
        inverse problems
        - Boundary conditions must be properly specified through
        the boundary dictionaries
    """

    def __init__(
        self,
        mesh,
        cells,
        boundary_points,
        cell_type: str,
        fe_order: int,
        fe_type: str,
        quad_order: int,
        quad_type: str,
        fe_transformation_type: str,
        bound_function_dict: dict,
        bound_condition_dict: dict,
        forcing_function,
        output_path: str,
    ) -> None:
        """
        The constructor of the Fespace2D class.

        Args:
            mesh: The mesh object.
            cells: The cells of the mesh.
            boundary_points: The boundary points of the mesh.
            cell_type: The type of the cell.
            fe_order: The order of the finite element.
            fe_type: The type of the finite element.
            quad_order: The order of the quadrature.
            quad_type: The type of the quadrature.
            fe_transformation_type: The type of the finite element transformation.
            bound_function_dict: The dictionary of the boundary functions.
            bound_condition_dict: The dictionary of the boundary conditions.
            forcing_function: The forcing function.
            output_path: The path to the output directory.

        Returns:
            None
        """
        self.mesh = mesh
        self.boundary_points = boundary_points
        self.cells = cells
        self.cell_type = cell_type
        self.fe_order = fe_order
        self.fe_type = fe_type
        self.quad_order = quad_order
        self.quad_type = quad_type

        self.fe_transformation_type = fe_transformation_type
        self.output_path = output_path
        self.bound_function_dict = bound_function_dict
        self.bound_condition_dict = bound_condition_dict
        self.forcing_function = forcing_function

    @abstractmethod
    def set_finite_elements(self) -> None:
        """
        Assigns the finite elements to each cell.

        This method initializes the finite element objects for each cell in the mesh.
        It creates an instance of the `FE2D_Cell` class for each cell, passing the necessary parameters.
        The finite element objects store information about the basis functions, gradients, Jacobians,
        quadrature points, weights, actual coordinates, and forcing functions associated with each cell.

        After initializing the finite element objects, this method prints the shape details of various matrices
        and updates the total number of degrees of freedom (dofs) for the entire mesh.

        Args:
            None

        Returns:
            None
        """

    @abstractmethod
    def generate_dirichlet_boundary_data(self) -> np.ndarray:
        """
        Generate Dirichlet boundary data.

        This function returns the boundary points and their corresponding values.

        Args:
            None

        Returns:
            np.ndarray: The boundary points and their values.

        Notes:
            The boundary points and values are stored in the `boundary_points` attribute of the `Fespace` object.
        """

    @abstractmethod
    def get_shape_function_val(self, cell_index: int) -> np.ndarray:
        """
        Get the actual values of the shape functions on a given cell.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: An array containing the actual values of the shape functions.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """

    @abstractmethod
    def get_shape_function_grad_x(self, cell_index: int) -> np.ndarray:
        """
        Get the gradient of the shape function with respect to the x-coordinate.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: An array containing the gradient of the shape function with respect to the x-coordinate.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """

    @abstractmethod
    def get_shape_function_grad_x_ref(self, cell_index: int) -> np.ndarray:
        """
        Get the gradient of the shape function with respect to the x-coordinate on the reference element.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: An array containing the gradient of the shape function with respect to the x-coordinate on the reference element.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """

    @abstractmethod
    def get_shape_function_grad_y(self, cell_index: int) -> np.ndarray:
        """
        Get the gradient of the shape function with respect to y at the given cell index.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: An array containing the gradient of the shape function with respect to y.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """

    @abstractmethod
    def get_shape_function_grad_y_ref(self, cell_index: int):
        """
        Get the gradient of the shape function with respect to y at the reference element.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: An array containing the gradient of the shape function with respect to y at the reference element.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.

        Notes:
            This function returns the gradient of the shape function with respect to y at the reference element
            for a given cell. The shape function gradient values are stored in the `basis_grady_at_quad_ref` array
            of the corresponding finite element cell. The `cell_index` parameter specifies the index of the cell
            for which the shape function gradient is required. If the `cell_index` is greater than the total number
            of cells, a `ValueError` is raised. The returned gradient values are copied from the `basis_grady_at_quad_ref` array to ensure immutability.
        """

    @abstractmethod
    def get_quadrature_actual_coordinates(self, cell_index: int) -> np.ndarray:
        """
        Get the actual coordinates of the quadrature points for a given cell.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: An array containing the actual coordinates of the quadrature points.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """

    @abstractmethod
    def get_forcing_function_values(self, cell_index: int) -> np.ndarray:
        """
        Get the forcing function values at the quadrature points.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: An array containing the forcing function values at the quadrature points.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.

        Notes:
            This function computes the forcing function values at the quadrature points for a given cell.
            It loops over all the basis functions and computes the integral using the actual coordinates
            and the basis functions at the quadrature points. The resulting values are stored in the
            `forcing_at_quad` attribute of the corresponding `fe_cell` object. The forcing function is evaluated using the `forcing_function` method of the `fe_cell`
            object.
        """

    @abstractmethod
    def get_sensor_data(self, exact_solution, num_points: int) -> np.ndarray:
        """
        Obtain sensor data (actual solution) at random points.

        Args:
            exact_solution (ndarray): The exact solution values.
            num_points (int): The number of points to sample from the domain.

        Returns:
            np.ndarray: The sensor data at the given points.

        Notes:
            This method is used in the inverse problem to obtain the sensor data at random points within the domain. Currently, it only works for problems with an analytical solution.
            Methodologies to obtain sensor data for problems from a file are not implemented yet.
            It is also not implemented for external or complex meshes.
        """

    @abstractmethod
    def get_sensor_data_external(
        self, exact_sol, num_points: int, file_name: str
    ) -> np.ndarray:
        """
        This method is used to obtain the sensor data from an external file when there is no analytical solution available.

        Args:
            exact_sol: The exact solution values.
            num_points: The number of points to sample from the domain.
            file_name: The name of the file containing the sensor data.

        Returns:
            np.ndarray: The sensor data at the given points based on the external file.

        """

__init__(mesh, cells, boundary_points, cell_type, fe_order, fe_type, quad_order, quad_type, fe_transformation_type, bound_function_dict, bound_condition_dict, forcing_function, output_path)

The constructor of the Fespace2D class.

Parameters:

Name Type Description Default
mesh

The mesh object.

 required
cells

The cells of the mesh.

 required
boundary_points

The boundary points of the mesh.

 required
cell_type str

The type of the cell.

 required
fe_order int

The order of the finite element.

 required
fe_type str

The type of the finite element.

 required
quad_order int

The order of the quadrature.

 required
quad_type str

The type of the quadrature.

 required
fe_transformation_type str

The type of the finite element transformation.

 required
bound_function_dict dict

The dictionary of the boundary functions.

 required
bound_condition_dict dict

The dictionary of the boundary conditions.

 required
forcing_function

The forcing function.

 required
output_path str

The path to the output directory.

 required

Returns:

Type Description
None

None
Source code in scirex/core/sciml/fe/fespace.py 
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
def __init__(
    self,
    mesh,
    cells,
    boundary_points,
    cell_type: str,
    fe_order: int,
    fe_type: str,
    quad_order: int,
    quad_type: str,
    fe_transformation_type: str,
    bound_function_dict: dict,
    bound_condition_dict: dict,
    forcing_function,
    output_path: str,
) -> None:
    """
    The constructor of the Fespace2D class.

    Args:
        mesh: The mesh object.
        cells: The cells of the mesh.
        boundary_points: The boundary points of the mesh.
        cell_type: The type of the cell.
        fe_order: The order of the finite element.
        fe_type: The type of the finite element.
        quad_order: The order of the quadrature.
        quad_type: The type of the quadrature.
        fe_transformation_type: The type of the finite element transformation.
        bound_function_dict: The dictionary of the boundary functions.
        bound_condition_dict: The dictionary of the boundary conditions.
        forcing_function: The forcing function.
        output_path: The path to the output directory.

    Returns:
        None
    """
    self.mesh = mesh
    self.boundary_points = boundary_points
    self.cells = cells
    self.cell_type = cell_type
    self.fe_order = fe_order
    self.fe_type = fe_type
    self.quad_order = quad_order
    self.quad_type = quad_type

    self.fe_transformation_type = fe_transformation_type
    self.output_path = output_path
    self.bound_function_dict = bound_function_dict
    self.bound_condition_dict = bound_condition_dict
    self.forcing_function = forcing_function

generate_dirichlet_boundary_data() abstractmethod

Generate Dirichlet boundary data.

This function returns the boundary points and their corresponding values.

Returns:

Type Description
ndarray

np.ndarray: The boundary points and their values.
Notes

The boundary points and values are stored in the boundary_points attribute of the Fespace object.

Source code in scirex/core/sciml/fe/fespace.py 
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
@abstractmethod
def generate_dirichlet_boundary_data(self) -> np.ndarray:
    """
    Generate Dirichlet boundary data.

    This function returns the boundary points and their corresponding values.

    Args:
        None

    Returns:
        np.ndarray: The boundary points and their values.

    Notes:
        The boundary points and values are stored in the `boundary_points` attribute of the `Fespace` object.
    """

get_forcing_function_values(cell_index) abstractmethod

Get the forcing function values at the quadrature points.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: An array containing the forcing function values at the quadrature points.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Notes

This function computes the forcing function values at the quadrature points for a given cell. It loops over all the basis functions and computes the integral using the actual coordinates and the basis functions at the quadrature points. The resulting values are stored in the forcing_at_quad attribute of the corresponding fe_cell object. The forcing function is evaluated using the forcing_function method of the fe_cell object.

Source code in scirex/core/sciml/fe/fespace.py 
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
@abstractmethod
def get_forcing_function_values(self, cell_index: int) -> np.ndarray:
    """
    Get the forcing function values at the quadrature points.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: An array containing the forcing function values at the quadrature points.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.

    Notes:
        This function computes the forcing function values at the quadrature points for a given cell.
        It loops over all the basis functions and computes the integral using the actual coordinates
        and the basis functions at the quadrature points. The resulting values are stored in the
        `forcing_at_quad` attribute of the corresponding `fe_cell` object. The forcing function is evaluated using the `forcing_function` method of the `fe_cell`
        object.
    """

get_quadrature_actual_coordinates(cell_index) abstractmethod

Get the actual coordinates of the quadrature points for a given cell.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: An array containing the actual coordinates of the quadrature points.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace.py 
280
281
282
283
284
285
286
287
288
289
290
291
292
293
@abstractmethod
def get_quadrature_actual_coordinates(self, cell_index: int) -> np.ndarray:
    """
    Get the actual coordinates of the quadrature points for a given cell.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: An array containing the actual coordinates of the quadrature points.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """

get_sensor_data(exact_solution, num_points) abstractmethod

Obtain sensor data (actual solution) at random points.

Parameters:

Name Type Description Default
exact_solution ndarray

The exact solution values.

 required
num_points int

The number of points to sample from the domain.

 required

Returns:

Type Description
ndarray

np.ndarray: The sensor data at the given points.
Notes

This method is used in the inverse problem to obtain the sensor data at random points within the domain. Currently, it only works for problems with an analytical solution. Methodologies to obtain sensor data for problems from a file are not implemented yet. It is also not implemented for external or complex meshes.

Source code in scirex/core/sciml/fe/fespace.py 
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
@abstractmethod
def get_sensor_data(self, exact_solution, num_points: int) -> np.ndarray:
    """
    Obtain sensor data (actual solution) at random points.

    Args:
        exact_solution (ndarray): The exact solution values.
        num_points (int): The number of points to sample from the domain.

    Returns:
        np.ndarray: The sensor data at the given points.

    Notes:
        This method is used in the inverse problem to obtain the sensor data at random points within the domain. Currently, it only works for problems with an analytical solution.
        Methodologies to obtain sensor data for problems from a file are not implemented yet.
        It is also not implemented for external or complex meshes.
    """

get_sensor_data_external(exact_sol, num_points, file_name) abstractmethod

This method is used to obtain the sensor data from an external file when there is no analytical solution available.

Parameters:

Name Type Description Default
exact_sol

The exact solution values.

 required
num_points int

The number of points to sample from the domain.

 required
file_name str

The name of the file containing the sensor data.

 required

Returns:

Type Description
ndarray

np.ndarray: The sensor data at the given points based on the external file.
Source code in scirex/core/sciml/fe/fespace.py 
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
@abstractmethod
def get_sensor_data_external(
    self, exact_sol, num_points: int, file_name: str
) -> np.ndarray:
    """
    This method is used to obtain the sensor data from an external file when there is no analytical solution available.

    Args:
        exact_sol: The exact solution values.
        num_points: The number of points to sample from the domain.
        file_name: The name of the file containing the sensor data.

    Returns:
        np.ndarray: The sensor data at the given points based on the external file.

    """

get_shape_function_grad_x(cell_index) abstractmethod

Get the gradient of the shape function with respect to the x-coordinate.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: An array containing the gradient of the shape function with respect to the x-coordinate.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace.py 
213
214
215
216
217
218
219
220
221
222
223
224
225
226
@abstractmethod
def get_shape_function_grad_x(self, cell_index: int) -> np.ndarray:
    """
    Get the gradient of the shape function with respect to the x-coordinate.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: An array containing the gradient of the shape function with respect to the x-coordinate.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """

get_shape_function_grad_x_ref(cell_index) abstractmethod

Get the gradient of the shape function with respect to the x-coordinate on the reference element.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: An array containing the gradient of the shape function with respect to the x-coordinate on the reference element.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace.py 
228
229
230
231
232
233
234
235
236
237
238
239
240
241
@abstractmethod
def get_shape_function_grad_x_ref(self, cell_index: int) -> np.ndarray:
    """
    Get the gradient of the shape function with respect to the x-coordinate on the reference element.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: An array containing the gradient of the shape function with respect to the x-coordinate on the reference element.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """

get_shape_function_grad_y(cell_index) abstractmethod

Get the gradient of the shape function with respect to y at the given cell index.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: An array containing the gradient of the shape function with respect to y.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace.py 
243
244
245
246
247
248
249
250
251
252
253
254
255
256
@abstractmethod
def get_shape_function_grad_y(self, cell_index: int) -> np.ndarray:
    """
    Get the gradient of the shape function with respect to y at the given cell index.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: An array containing the gradient of the shape function with respect to y.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """

get_shape_function_grad_y_ref(cell_index) abstractmethod

Get the gradient of the shape function with respect to y at the reference element.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description

np.ndarray: An array containing the gradient of the shape function with respect to y at the reference element.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Notes

This function returns the gradient of the shape function with respect to y at the reference element for a given cell. The shape function gradient values are stored in the basis_grady_at_quad_ref array of the corresponding finite element cell. The cell_index parameter specifies the index of the cell for which the shape function gradient is required. If the cell_index is greater than the total number of cells, a ValueError is raised. The returned gradient values are copied from the basis_grady_at_quad_ref array to ensure immutability.

Source code in scirex/core/sciml/fe/fespace.py 
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
@abstractmethod
def get_shape_function_grad_y_ref(self, cell_index: int):
    """
    Get the gradient of the shape function with respect to y at the reference element.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: An array containing the gradient of the shape function with respect to y at the reference element.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.

    Notes:
        This function returns the gradient of the shape function with respect to y at the reference element
        for a given cell. The shape function gradient values are stored in the `basis_grady_at_quad_ref` array
        of the corresponding finite element cell. The `cell_index` parameter specifies the index of the cell
        for which the shape function gradient is required. If the `cell_index` is greater than the total number
        of cells, a `ValueError` is raised. The returned gradient values are copied from the `basis_grady_at_quad_ref` array to ensure immutability.
    """

get_shape_function_val(cell_index) abstractmethod

Get the actual values of the shape functions on a given cell.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: An array containing the actual values of the shape functions.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace.py 
198
199
200
201
202
203
204
205
206
207
208
209
210
211
@abstractmethod
def get_shape_function_val(self, cell_index: int) -> np.ndarray:
    """
    Get the actual values of the shape functions on a given cell.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: An array containing the actual values of the shape functions.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """

set_finite_elements() abstractmethod

Assigns the finite elements to each cell.

This method initializes the finite element objects for each cell in the mesh. It creates an instance of the FE2D_Cell class for each cell, passing the necessary parameters. The finite element objects store information about the basis functions, gradients, Jacobians, quadrature points, weights, actual coordinates, and forcing functions associated with each cell.

After initializing the finite element objects, this method prints the shape details of various matrices and updates the total number of degrees of freedom (dofs) for the entire mesh.

Returns:

Type Description
None

None
Source code in scirex/core/sciml/fe/fespace.py 
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
@abstractmethod
def set_finite_elements(self) -> None:
    """
    Assigns the finite elements to each cell.

    This method initializes the finite element objects for each cell in the mesh.
    It creates an instance of the `FE2D_Cell` class for each cell, passing the necessary parameters.
    The finite element objects store information about the basis functions, gradients, Jacobians,
    quadrature points, weights, actual coordinates, and forcing functions associated with each cell.

    After initializing the finite element objects, this method prints the shape details of various matrices
    and updates the total number of degrees of freedom (dofs) for the entire mesh.

    Args:
        None

    Returns:
        None
    """

Fespace2D

Bases: Fespace

Represents a finite element space in 2D. This class provides functionality for handling 2D finite element spaces, including mesh generation, basis function evaluation, and boundary condition handling.

Parameters:

Name Type Description Default
mesh Mesh

The mesh object containing the mesh information.

 required
cells ndarray

The cell information from the mesh.

 required
boundary_points dict

The boundary points information from the mesh.

 required
cell_type str

The type of the cell (e.g., 'quadrilateral').

 required
fe_order int

The order of the finite element basis functions.

 required
fe_type str

The type of the finite element basis functions (e.g., 'legendre').

 required
quad_order int

The order of the quadrature rule.

 required
quad_type str

The type of the quadrature rule (e.g., 'gauss-legendre').

 required
fe_transformation_type str

The type of the finite element transformation (e.g., 'affine').

 required
bound_function_dict dict

A dictionary containing the boundary functions.

 required
bound_condition_dict dict

A dictionary containing the boundary conditions.

 required
forcing_function function

The forcing function for the problem.

 required
output_path str

The path to save the output files.

 required
generate_mesh_plot bool

Flag to generate the mesh plot (default: False).

 False

Raises:

Type Description
ValueError

If the cell type is not supported.

Returns:

Type Description

None
Source code in scirex/core/sciml/fe/fespace2d.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
class Fespace2D(Fespace):
    """
    Represents a finite element space in 2D. This class provides functionality for handling 2D finite element spaces,
    including mesh generation, basis function evaluation, and boundary condition handling.

    Args:
        mesh (meshio.Mesh): The mesh object containing the mesh information.
        cells (np.ndarray): The cell information from the mesh.
        boundary_points (dict): The boundary points information from the mesh.
        cell_type (str): The type of the cell (e.g., 'quadrilateral').
        fe_order (int): The order of the finite element basis functions.
        fe_type (str): The type of the finite element basis functions (e.g., 'legendre').
        quad_order (int): The order of the quadrature rule.
        quad_type (str): The type of the quadrature rule (e.g., 'gauss-legendre').
        fe_transformation_type (str): The type of the finite element transformation (e.g., 'affine').
        bound_function_dict (dict): A dictionary containing the boundary functions.
        bound_condition_dict (dict): A dictionary containing the boundary conditions.
        forcing_function (function): The forcing function for the problem.
        output_path (str): The path to save the output files.
        generate_mesh_plot (bool): Flag to generate the mesh plot (default: False).

    Raises:
        ValueError: If the cell type is not supported.

    Returns:
        None
    """

    def __init__(
        self,
        mesh,
        cells,
        boundary_points,
        cell_type: str,
        fe_order: int,
        fe_type: str,
        quad_order: int,
        quad_type: str,
        fe_transformation_type: str,
        bound_function_dict: dict,
        bound_condition_dict: dict,
        forcing_function,
        output_path: str,
        generate_mesh_plot: bool = False,
    ) -> None:
        """
        The constructor of the Fespace2D class.
        """
        # call the constructor of the parent class
        super().__init__(
            mesh=mesh,
            cells=cells,
            boundary_points=boundary_points,
            cell_type=cell_type,
            fe_order=fe_order,
            fe_type=fe_type,
            quad_order=quad_order,
            quad_type=quad_type,
            fe_transformation_type=fe_transformation_type,
            bound_function_dict=bound_function_dict,
            bound_condition_dict=bound_condition_dict,
            forcing_function=forcing_function,
            output_path=output_path,
        )

        if self.cell_type == "triangle":
            raise ValueError(
                "Triangle Mesh is not supported yet"
            )  # added by thivin - to remove support for triangular mesh

        self.generate_mesh_plot = generate_mesh_plot

        # to be calculated in the plot function
        self.total_dofs = 0
        self.total_boundary_dofs = 0

        # to be calculated on get_boundary_data_dirichlet function
        self.total_dirichlet_dofs = 0

        # get the number of cells
        self.n_cells = self.cells.shape[0]

        self.fe_cell = []

        # Function which assigns the fe_cell for each cell
        self.set_finite_elements()

        # generate the plot of the mesh
        if self.generate_mesh_plot:
            self.generate_plot(self.output_path)
        # self.generate_plot(self.output_path)

        # Obtain boundary Data
        self.dirichlet_boundary_data = self.generate_dirichlet_boundary_data()

        title = [
            "Number of Cells",
            "Number of Quadrature Points",
            "Number of Dirichlet Boundary Points",
            "Quadrature Order",
            "fe Order",
            "fe Type",
            "fe Transformation Type",
        ]
        values = [
            self.n_cells,
            self.total_dofs,
            self.total_dirichlet_dofs,
            self.quad_order,
            self.fe_order,
            self.fe_type,
            self.fe_transformation_type,
        ]
        # print the table
        print_table("fe Space Information", ["Property", "Value"], title, values)

    def set_finite_elements(self) -> None:
        """
        Assigns the finite elements to each cell.

        This method initializes the finite element objects for each cell in the mesh.
        It creates an instance of the `FE2D_Cell` class for each cell, passing the necessary parameters.
        The finite element objects store information about the basis functions, gradients, Jacobians,
        quadrature points, weights, actual coordinates, and forcing functions associated with each cell.

        After initializing the finite element objects, this method prints the shape details of various matrices
        and updates the total number of degrees of freedom (dofs) for the entire mesh.

        :return: None
        """
        progress_bar = tqdm(
            total=self.n_cells,
            desc="Fe2D_cell Setup",
            unit="cells_assembled",
            bar_format="{l_bar}{bar:40}{r_bar}{bar:-10b}",
            colour="blue",
            ncols=100,
        )

        dof = 0
        for i in range(self.n_cells):
            self.fe_cell.append(
                FE2D_Cell(
                    self.cells[i],
                    self.cell_type,
                    self.fe_order,
                    self.fe_type,
                    self.quad_order,
                    self.quad_type,
                    self.fe_transformation_type,
                    self.forcing_function,
                )
            )

            # obtain the shape of the basis function (n_test, N_quad)
            dof += self.fe_cell[i].basis_at_quad.shape[1]

            progress_bar.update(1)
        # print the Shape details of all the matrices from cell 0 using print_table function
        title = [
            "Shape function Matrix Shape",
            "Shape function Gradient Matrix Shape",
            "Jacobian Matrix Shape",
            "Quadrature Points Shape",
            "Quadrature Weights Shape",
            "Quadrature Actual Coordinates Shape",
            "Forcing Function Shape",
        ]
        values = [
            self.fe_cell[0].basis_at_quad.shape,
            self.fe_cell[0].basis_gradx_at_quad.shape,
            self.fe_cell[0].jacobian.shape,
            self.fe_cell[0].quad_xi.shape,
            self.fe_cell[0].quad_weight.shape,
            self.fe_cell[0].quad_actual_coordinates.shape,
            self.fe_cell[0].forcing_at_quad.shape,
        ]
        print_table("fe Matrix Shapes", ["Matrix", "Shape"], title, values)

        # update the total number of dofs
        self.total_dofs = dof

    def generate_plot(self, output_path) -> None:
        """
        Generate a plot of the mesh.

        Args:
            output_path (str): The path to save the output files.

        Returns:
            None
        """
        total_quad = 0
        marker_list = [
            "o",
            ".",
            ",",
            "x",
            "+",
            "P",
            "s",
            "D",
            "d",
            "^",
            "v",
            "<",
            ">",
            "p",
            "h",
            "H",
        ]

        print(f"[INFO] : Generating the plot of the mesh")
        # Plot the mesh
        plt.figure(figsize=(6.4, 4.8), dpi=300)

        # label flag ( to add the label only once)
        label_set = False

        # plot every cell as a quadrilateral
        # loop over all the cells
        for i in range(self.n_cells):
            # get the coordinates of the cell
            x = self.fe_cell[i].cell_coordinates[:, 0]
            y = self.fe_cell[i].cell_coordinates[:, 1]

            # add the first point to the end of the array
            x = np.append(x, x[0])
            y = np.append(y, y[0])

            plt.plot(x, y, "k-", linewidth=0.5)

            # plot the quadrature points
            x_quad = self.fe_cell[i].quad_actual_coordinates[:, 0]
            y_quad = self.fe_cell[i].quad_actual_coordinates[:, 1]

            total_quad += x_quad.shape[0]

            if not label_set:
                plt.scatter(
                    x_quad, y_quad, marker="x", color="b", s=2, label="Quad Pts"
                )
                label_set = True
            else:
                plt.scatter(x_quad, y_quad, marker="x", color="b", s=2)

        self.total_dofs = total_quad

        bound_dof = 0
        # plot the boundary points
        # loop over all the boundary tags
        for i, (bound_id, bound_pts) in enumerate(self.boundary_points.items()):
            # get the coordinates of the boundary points
            x = bound_pts[:, 0]
            y = bound_pts[:, 1]

            # add the first point to the end of the array
            x = np.append(x, x[0])
            y = np.append(y, y[0])

            bound_dof += x.shape[0]

            plt.scatter(
                x, y, marker=marker_list[i + 1], s=2, label=f"Bd-id : {bound_id}"
            )

        self.total_boundary_dofs = bound_dof

        plt.legend(bbox_to_anchor=(0.85, 1.02))
        plt.axis("equal")
        plt.axis("off")
        plt.tight_layout()

        plt.savefig(str(Path(output_path) / "mesh.png"), bbox_inches="tight")
        plt.savefig(str(Path(output_path) / "mesh.svg"), bbox_inches="tight")

        # print the total number of quadrature points
        print(f"Plots generated")
        print(f"[INFO] : Total number of cells = {self.n_cells}")
        print(f"[INFO] : Total number of quadrature points = {self.total_dofs}")
        print(f"[INFO] : Total number of boundary points = {self.total_boundary_dofs}")

    def generate_dirichlet_boundary_data(self) -> np.ndarray:
        """
        Generate Dirichlet boundary data. This function returns the boundary points and their corresponding values.

        Args:
            None

        Returns:
            tuple: The boundary points and their values as numpy arrays.
        """
        x = []
        y = []
        for bound_id, bound_pts in self.boundary_points.items():
            # get the coordinates of the boundary points
            for pt in bound_pts:
                pt_new = np.array([pt[0], pt[1]], dtype=np.float64)
                x.append(pt_new)
                val = np.array(
                    self.bound_function_dict[bound_id](pt[0], pt[1]), dtype=np.float64
                ).reshape(-1, 1)
                y.append(val)

        print(f"[INFO] : Total number of Dirichlet boundary points = {len(x)}")
        self.total_dirichlet_dofs = len(x)
        print(f"[INFO] : Shape of Dirichlet-X = {np.array(x).shape}")
        print(f"[INFO] : Shape of Y = {np.array(y).shape}")

        return x, y

    def generate_dirichlet_boundary_data_vector(self, component: int) -> np.ndarray:
        """
        Generate the boundary data vector for the Dirichlet boundary condition. This function returns the boundary points and their corresponding values for a specific component.

        Args:
            component (int): The component of the boundary data vector.

        Returns:
            tuple: The boundary points and their values as numpy arrays.
        """
        x = []
        y = []
        for bound_id, bound_pts in self.boundary_points.items():
            # get the coordinates of the boundary points
            for pt in bound_pts:
                pt_new = np.array([pt[0], pt[1]], dtype=np.float64)
                x.append(pt_new)
                val = np.array(
                    self.bound_function_dict[bound_id](pt[0], pt[1])[component],
                    dtype=np.float64,
                ).reshape(-1, 1)
                y.append(val)

        return x, y

    def get_shape_function_val(self, cell_index: int) -> np.ndarray:
        """
        Get the actual values of the shape functions on a given cell.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: The actual values of the shape functions on the given cell.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """
        if cell_index >= len(self.fe_cell) or cell_index < 0:
            raise ValueError(
                f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
            )

        return self.fe_cell[cell_index].basis_at_quad.copy()

    def get_shape_function_grad_x(self, cell_index: int) -> np.ndarray:
        """
        Get the gradient of the shape function with respect to the x-coordinate.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: The actual values of the shape functions on the given cell.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """
        if cell_index >= len(self.fe_cell) or cell_index < 0:
            raise ValueError(
                f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
            )

        return self.fe_cell[cell_index].basis_gradx_at_quad.copy()

    def get_shape_function_grad_x_ref(self, cell_index: int) -> np.ndarray:
        """
        Get the gradient of the shape function with respect to the x-coordinate on the reference element.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: The actual values of the shape functions on the given cell.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """
        if cell_index >= len(self.fe_cell) or cell_index < 0:
            raise ValueError(
                f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
            )

        return self.fe_cell[cell_index].basis_gradx_at_quad_ref.copy()

    def get_shape_function_grad_y(self, cell_index: int) -> np.ndarray:
        """
        Get the gradient of the shape function with respect to y at the given cell index.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: The actual values of the shape functions on the given cell.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """
        if cell_index >= len(self.fe_cell) or cell_index < 0:
            raise ValueError(
                f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
            )

        return self.fe_cell[cell_index].basis_grady_at_quad.copy()

    def get_shape_function_grad_y_ref(self, cell_index: int):
        """
        Get the gradient of the shape function with respect to y at the reference element.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: The actual values of the shape functions on the given cell.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.

        Note:
            This function returns the gradient of the shape function with respect to y at the reference element
            for a given cell. The shape function gradient values are stored in the `basis_grady_at_quad_ref` array
            of the corresponding finite element cell. The `cell_index` parameter specifies the index of the cell
            for which the shape function gradient is required. If the `cell_index` is greater than the total number
            of cells, a `ValueError` is raised. The returned gradient values are copied from the `basis_grady_at_quad_ref` array to ensure immutability.
        """
        if cell_index >= len(self.fe_cell) or cell_index < 0:
            raise ValueError(
                f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
            )

        return self.fe_cell[cell_index].basis_grady_at_quad_ref.copy()

    def get_quadrature_actual_coordinates(self, cell_index: int) -> np.ndarray:
        """
        Get the actual coordinates of the quadrature points for a given cell.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: An array containing the actual coordinates of the quadrature points.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """
        if cell_index >= len(self.fe_cell) or cell_index < 0:
            raise ValueError(
                f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
            )

        return self.fe_cell[cell_index].quad_actual_coordinates.copy()

    def get_quadrature_weights(self, cell_index: int) -> np.ndarray:
        """
        Return the quadrature weights for a given cell.

        Args:
            cell_index (int): The index of the cell for which the quadrature weights are needed.

        Returns:
            np.ndarray: The quadrature weights for the given cell  of dimension (N_Quad_Points, 1).

        Raises:
            ValueError: If cell_index is greater than the number of cells.
        """
        if cell_index >= len(self.fe_cell) or cell_index < 0:
            raise ValueError(
                f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
            )

        return self.fe_cell[cell_index].mult.copy()

    def get_forcing_function_values(self, cell_index: int) -> np.ndarray:
        """
        Get the forcing function values at the quadrature points.

        Args:
            cell_index (int): The index of the cell.

        Returns:
            np.ndarray: The forcing function values at the quadrature points.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.

        Note:
            This function computes the forcing function values at the quadrature points for a given cell.
            It loops over all the basis functions and computes the integral using the actual coordinates
            and the basis functions at the quadrature points. The resulting values are stored in the
            `forcing_at_quad` attribute of the corresponding `fe_cell` object.
        """
        if cell_index >= len(self.fe_cell) or cell_index < 0:
            raise ValueError(
                f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
            )

        # Changed by Thivin: To assemble the forcing function at the quadrature points here in the fespace
        # so that it can be used to handle multiple dimensions on a vector valud problem

        # get number of shape functions
        n_shape_functions = self.fe_cell[cell_index].basis_function.num_shape_functions

        # Loop over all the basis functions and compute the integral
        f_integral = np.zeros((n_shape_functions, 1), dtype=np.float64)

        for i in range(n_shape_functions):
            val = 0
            for q in range(self.fe_cell[cell_index].basis_at_quad.shape[1]):
                x = self.fe_cell[cell_index].quad_actual_coordinates[q, 0]
                y = self.fe_cell[cell_index].quad_actual_coordinates[q, 1]
                # print("f_values[q] = ",f_values[q])

                # the Jacobian and the quadrature weights are pre multiplied to the basis functions
                val += (self.fe_cell[cell_index].basis_at_quad[i, q]) * self.fe_cell[
                    cell_index
                ].forcing_function(x, y)
                # print("val = ", val)

            f_integral[i] = val

        self.fe_cell[cell_index].forcing_at_quad = f_integral

        return self.fe_cell[cell_index].forcing_at_quad.copy()

    def get_forcing_function_values_vector(
        self, cell_index: int, component: int
    ) -> np.ndarray:
        """
        This function will return the forcing function values at the quadrature points
        based on the Component of the RHS Needed, for vector valued problems

        Args:
            cell_index (int): The index of the cell.
            component (int): The component of the forcing function.

        Returns:
            np.ndarray: The forcing function values at the quadrature points.

        Raises:
            ValueError: If the cell_index is greater than the number of cells.
        """
        if cell_index >= len(self.fe_cell) or cell_index < 0:
            raise ValueError(
                f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
            )

        # get the coordinates
        x = self.fe_cell[cell_index].quad_actual_coordinates[:, 0]
        y = self.fe_cell[cell_index].quad_actual_coordinates[:, 1]

        # compute the forcing function values
        f_values = self.fe_cell[cell_index].forcing_function(x, y)[component]

        # compute the integral
        f_integral = np.sum(self.fe_cell[cell_index].basis_at_quad * f_values, axis=1)

        self.fe_cell[cell_index].forcing_at_quad = f_integral.reshape(-1, 1)

        return self.fe_cell[cell_index].forcing_at_quad.copy()

    def get_sensor_data(self, exact_solution, num_points: int):
        """
        Obtain sensor data (actual solution) at random points.

        This method is used in the inverse problem to obtain the sensor data at random points within the domain.
        Currently, it only works for problems with an analytical solution.
        Methodologies to obtain sensor data for problems from a file are not implemented yet.
        It is also not implemented for external or complex meshes.

        Args:
            exact_solution (function): The exact solution function.
            num_points (int): The number of points to sample.

        Returns:
            Tuple: A tuple containing two arrays: sensor points and the exact solution values.
        """
        # generate random points within the bounds of the domain
        # get the bounds of the domain
        x_min = np.min(self.mesh.points[:, 0])
        x_max = np.max(self.mesh.points[:, 0])
        y_min = np.min(self.mesh.points[:, 1])
        y_max = np.max(self.mesh.points[:, 1])
        # sample n random points within the bounds of the domain
        # Generate points in the unit square

        num_internal_points = int(num_points * 0.9)

        points = lhs(2, samples=num_internal_points)
        points[:, 0] = x_min + (x_max - x_min) * points[:, 0]
        points[:, 1] = y_min + (y_max - y_min) * points[:, 1]
        # get the exact solution at the points
        exact_sol = exact_solution(points[:, 0], points[:, 1])

        # print the shape of the points and the exact solution
        print(f"[INFO] : Number of sensor points = {points.shape[0]}")
        print(f"[INFO] : Shape of sensor points = {points.shape}")

        # plot the points
        plt.figure(figsize=(6.4, 4.8), dpi=300)
        plt.scatter(points[:, 0], points[:, 1], marker="x", color="r", s=2)
        plt.axis("equal")
        plt.title("Sensor Points")
        plt.tight_layout()
        plt.savefig("sensor_points.png", bbox_inches="tight")

        return points, exact_sol

    def get_sensor_data_external(self, exact_sol, num_points: int, file_name: str):
        """
        This method is used to obtain the sensor data from an external file.

        Args:
            exact_sol (function): The exact solution function.
            num_points (int): The number of points to sample.
            file_name (str): The name of the file containing the sensor data.

        Returns:
            Tuple: A tuple containing two arrays: sensor points and the exact solution values.

        Note:
            This method reads the sensor data from a file and samples `num_points` from the data.
            The sensor data is then returned as a tuple containing the sensor points and the exact solution values.
        """
        # use pandas to read the file
        df = pd.read_csv(file_name)

        x = df.iloc[:, 0].values
        y = df.iloc[:, 1].values
        exact_sol = df.iloc[:, 2].values

        # now sample num_points from the data
        indices = np.random.randint(0, x.shape[0], num_points)

        x = x[indices]
        y = y[indices]
        exact_sol = exact_sol[indices]

        # stack them together
        points = np.stack((x, y), axis=1)

        return points, exact_sol

__init__(mesh, cells, boundary_points, cell_type, fe_order, fe_type, quad_order, quad_type, fe_transformation_type, bound_function_dict, bound_condition_dict, forcing_function, output_path, generate_mesh_plot=False)

The constructor of the Fespace2D class.

Source code in scirex/core/sciml/fe/fespace2d.py 
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
def __init__(
    self,
    mesh,
    cells,
    boundary_points,
    cell_type: str,
    fe_order: int,
    fe_type: str,
    quad_order: int,
    quad_type: str,
    fe_transformation_type: str,
    bound_function_dict: dict,
    bound_condition_dict: dict,
    forcing_function,
    output_path: str,
    generate_mesh_plot: bool = False,
) -> None:
    """
    The constructor of the Fespace2D class.
    """
    # call the constructor of the parent class
    super().__init__(
        mesh=mesh,
        cells=cells,
        boundary_points=boundary_points,
        cell_type=cell_type,
        fe_order=fe_order,
        fe_type=fe_type,
        quad_order=quad_order,
        quad_type=quad_type,
        fe_transformation_type=fe_transformation_type,
        bound_function_dict=bound_function_dict,
        bound_condition_dict=bound_condition_dict,
        forcing_function=forcing_function,
        output_path=output_path,
    )

    if self.cell_type == "triangle":
        raise ValueError(
            "Triangle Mesh is not supported yet"
        )  # added by thivin - to remove support for triangular mesh

    self.generate_mesh_plot = generate_mesh_plot

    # to be calculated in the plot function
    self.total_dofs = 0
    self.total_boundary_dofs = 0

    # to be calculated on get_boundary_data_dirichlet function
    self.total_dirichlet_dofs = 0

    # get the number of cells
    self.n_cells = self.cells.shape[0]

    self.fe_cell = []

    # Function which assigns the fe_cell for each cell
    self.set_finite_elements()

    # generate the plot of the mesh
    if self.generate_mesh_plot:
        self.generate_plot(self.output_path)
    # self.generate_plot(self.output_path)

    # Obtain boundary Data
    self.dirichlet_boundary_data = self.generate_dirichlet_boundary_data()

    title = [
        "Number of Cells",
        "Number of Quadrature Points",
        "Number of Dirichlet Boundary Points",
        "Quadrature Order",
        "fe Order",
        "fe Type",
        "fe Transformation Type",
    ]
    values = [
        self.n_cells,
        self.total_dofs,
        self.total_dirichlet_dofs,
        self.quad_order,
        self.fe_order,
        self.fe_type,
        self.fe_transformation_type,
    ]
    # print the table
    print_table("fe Space Information", ["Property", "Value"], title, values)

generate_dirichlet_boundary_data()

Generate Dirichlet boundary data. This function returns the boundary points and their corresponding values.

Returns:

Name Type Description
tuple ndarray

The boundary points and their values as numpy arrays.
Source code in scirex/core/sciml/fe/fespace2d.py 
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
def generate_dirichlet_boundary_data(self) -> np.ndarray:
    """
    Generate Dirichlet boundary data. This function returns the boundary points and their corresponding values.

    Args:
        None

    Returns:
        tuple: The boundary points and their values as numpy arrays.
    """
    x = []
    y = []
    for bound_id, bound_pts in self.boundary_points.items():
        # get the coordinates of the boundary points
        for pt in bound_pts:
            pt_new = np.array([pt[0], pt[1]], dtype=np.float64)
            x.append(pt_new)
            val = np.array(
                self.bound_function_dict[bound_id](pt[0], pt[1]), dtype=np.float64
            ).reshape(-1, 1)
            y.append(val)

    print(f"[INFO] : Total number of Dirichlet boundary points = {len(x)}")
    self.total_dirichlet_dofs = len(x)
    print(f"[INFO] : Shape of Dirichlet-X = {np.array(x).shape}")
    print(f"[INFO] : Shape of Y = {np.array(y).shape}")

    return x, y

generate_dirichlet_boundary_data_vector(component)

Generate the boundary data vector for the Dirichlet boundary condition. This function returns the boundary points and their corresponding values for a specific component.

Parameters:

Name Type Description Default
component int

The component of the boundary data vector.

 required

Returns:

Name Type Description
tuple ndarray

The boundary points and their values as numpy arrays.
Source code in scirex/core/sciml/fe/fespace2d.py 
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
def generate_dirichlet_boundary_data_vector(self, component: int) -> np.ndarray:
    """
    Generate the boundary data vector for the Dirichlet boundary condition. This function returns the boundary points and their corresponding values for a specific component.

    Args:
        component (int): The component of the boundary data vector.

    Returns:
        tuple: The boundary points and their values as numpy arrays.
    """
    x = []
    y = []
    for bound_id, bound_pts in self.boundary_points.items():
        # get the coordinates of the boundary points
        for pt in bound_pts:
            pt_new = np.array([pt[0], pt[1]], dtype=np.float64)
            x.append(pt_new)
            val = np.array(
                self.bound_function_dict[bound_id](pt[0], pt[1])[component],
                dtype=np.float64,
            ).reshape(-1, 1)
            y.append(val)

    return x, y

generate_plot(output_path)

Generate a plot of the mesh.

Parameters:

Name Type Description Default
output_path str

The path to save the output files.

 required

Returns:

Type Description
None

None
Source code in scirex/core/sciml/fe/fespace2d.py 
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
def generate_plot(self, output_path) -> None:
    """
    Generate a plot of the mesh.

    Args:
        output_path (str): The path to save the output files.

    Returns:
        None
    """
    total_quad = 0
    marker_list = [
        "o",
        ".",
        ",",
        "x",
        "+",
        "P",
        "s",
        "D",
        "d",
        "^",
        "v",
        "<",
        ">",
        "p",
        "h",
        "H",
    ]

    print(f"[INFO] : Generating the plot of the mesh")
    # Plot the mesh
    plt.figure(figsize=(6.4, 4.8), dpi=300)

    # label flag ( to add the label only once)
    label_set = False

    # plot every cell as a quadrilateral
    # loop over all the cells
    for i in range(self.n_cells):
        # get the coordinates of the cell
        x = self.fe_cell[i].cell_coordinates[:, 0]
        y = self.fe_cell[i].cell_coordinates[:, 1]

        # add the first point to the end of the array
        x = np.append(x, x[0])
        y = np.append(y, y[0])

        plt.plot(x, y, "k-", linewidth=0.5)

        # plot the quadrature points
        x_quad = self.fe_cell[i].quad_actual_coordinates[:, 0]
        y_quad = self.fe_cell[i].quad_actual_coordinates[:, 1]

        total_quad += x_quad.shape[0]

        if not label_set:
            plt.scatter(
                x_quad, y_quad, marker="x", color="b", s=2, label="Quad Pts"
            )
            label_set = True
        else:
            plt.scatter(x_quad, y_quad, marker="x", color="b", s=2)

    self.total_dofs = total_quad

    bound_dof = 0
    # plot the boundary points
    # loop over all the boundary tags
    for i, (bound_id, bound_pts) in enumerate(self.boundary_points.items()):
        # get the coordinates of the boundary points
        x = bound_pts[:, 0]
        y = bound_pts[:, 1]

        # add the first point to the end of the array
        x = np.append(x, x[0])
        y = np.append(y, y[0])

        bound_dof += x.shape[0]

        plt.scatter(
            x, y, marker=marker_list[i + 1], s=2, label=f"Bd-id : {bound_id}"
        )

    self.total_boundary_dofs = bound_dof

    plt.legend(bbox_to_anchor=(0.85, 1.02))
    plt.axis("equal")
    plt.axis("off")
    plt.tight_layout()

    plt.savefig(str(Path(output_path) / "mesh.png"), bbox_inches="tight")
    plt.savefig(str(Path(output_path) / "mesh.svg"), bbox_inches="tight")

    # print the total number of quadrature points
    print(f"Plots generated")
    print(f"[INFO] : Total number of cells = {self.n_cells}")
    print(f"[INFO] : Total number of quadrature points = {self.total_dofs}")
    print(f"[INFO] : Total number of boundary points = {self.total_boundary_dofs}")

get_forcing_function_values(cell_index)

Get the forcing function values at the quadrature points.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: The forcing function values at the quadrature points.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Note

This function computes the forcing function values at the quadrature points for a given cell. It loops over all the basis functions and computes the integral using the actual coordinates and the basis functions at the quadrature points. The resulting values are stored in the forcing_at_quad attribute of the corresponding fe_cell object.

Source code in scirex/core/sciml/fe/fespace2d.py 
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
def get_forcing_function_values(self, cell_index: int) -> np.ndarray:
    """
    Get the forcing function values at the quadrature points.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: The forcing function values at the quadrature points.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.

    Note:
        This function computes the forcing function values at the quadrature points for a given cell.
        It loops over all the basis functions and computes the integral using the actual coordinates
        and the basis functions at the quadrature points. The resulting values are stored in the
        `forcing_at_quad` attribute of the corresponding `fe_cell` object.
    """
    if cell_index >= len(self.fe_cell) or cell_index < 0:
        raise ValueError(
            f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
        )

    # Changed by Thivin: To assemble the forcing function at the quadrature points here in the fespace
    # so that it can be used to handle multiple dimensions on a vector valud problem

    # get number of shape functions
    n_shape_functions = self.fe_cell[cell_index].basis_function.num_shape_functions

    # Loop over all the basis functions and compute the integral
    f_integral = np.zeros((n_shape_functions, 1), dtype=np.float64)

    for i in range(n_shape_functions):
        val = 0
        for q in range(self.fe_cell[cell_index].basis_at_quad.shape[1]):
            x = self.fe_cell[cell_index].quad_actual_coordinates[q, 0]
            y = self.fe_cell[cell_index].quad_actual_coordinates[q, 1]
            # print("f_values[q] = ",f_values[q])

            # the Jacobian and the quadrature weights are pre multiplied to the basis functions
            val += (self.fe_cell[cell_index].basis_at_quad[i, q]) * self.fe_cell[
                cell_index
            ].forcing_function(x, y)
            # print("val = ", val)

        f_integral[i] = val

    self.fe_cell[cell_index].forcing_at_quad = f_integral

    return self.fe_cell[cell_index].forcing_at_quad.copy()

get_forcing_function_values_vector(cell_index, component)

This function will return the forcing function values at the quadrature points based on the Component of the RHS Needed, for vector valued problems

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required
component int

The component of the forcing function.

 required

Returns:

Type Description
ndarray

np.ndarray: The forcing function values at the quadrature points.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace2d.py 
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
def get_forcing_function_values_vector(
    self, cell_index: int, component: int
) -> np.ndarray:
    """
    This function will return the forcing function values at the quadrature points
    based on the Component of the RHS Needed, for vector valued problems

    Args:
        cell_index (int): The index of the cell.
        component (int): The component of the forcing function.

    Returns:
        np.ndarray: The forcing function values at the quadrature points.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """
    if cell_index >= len(self.fe_cell) or cell_index < 0:
        raise ValueError(
            f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
        )

    # get the coordinates
    x = self.fe_cell[cell_index].quad_actual_coordinates[:, 0]
    y = self.fe_cell[cell_index].quad_actual_coordinates[:, 1]

    # compute the forcing function values
    f_values = self.fe_cell[cell_index].forcing_function(x, y)[component]

    # compute the integral
    f_integral = np.sum(self.fe_cell[cell_index].basis_at_quad * f_values, axis=1)

    self.fe_cell[cell_index].forcing_at_quad = f_integral.reshape(-1, 1)

    return self.fe_cell[cell_index].forcing_at_quad.copy()

get_quadrature_actual_coordinates(cell_index)

Get the actual coordinates of the quadrature points for a given cell.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: An array containing the actual coordinates of the quadrature points.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace2d.py 
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
def get_quadrature_actual_coordinates(self, cell_index: int) -> np.ndarray:
    """
    Get the actual coordinates of the quadrature points for a given cell.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: An array containing the actual coordinates of the quadrature points.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """
    if cell_index >= len(self.fe_cell) or cell_index < 0:
        raise ValueError(
            f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
        )

    return self.fe_cell[cell_index].quad_actual_coordinates.copy()

get_quadrature_weights(cell_index)

Return the quadrature weights for a given cell.

Parameters:

Name Type Description Default
cell_index int

The index of the cell for which the quadrature weights are needed.

 required

Returns:

Type Description
ndarray

np.ndarray: The quadrature weights for the given cell of dimension (N_Quad_Points, 1).

Raises:

Type Description
ValueError

If cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace2d.py 
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
def get_quadrature_weights(self, cell_index: int) -> np.ndarray:
    """
    Return the quadrature weights for a given cell.

    Args:
        cell_index (int): The index of the cell for which the quadrature weights are needed.

    Returns:
        np.ndarray: The quadrature weights for the given cell  of dimension (N_Quad_Points, 1).

    Raises:
        ValueError: If cell_index is greater than the number of cells.
    """
    if cell_index >= len(self.fe_cell) or cell_index < 0:
        raise ValueError(
            f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
        )

    return self.fe_cell[cell_index].mult.copy()

get_sensor_data(exact_solution, num_points)

Obtain sensor data (actual solution) at random points.

This method is used in the inverse problem to obtain the sensor data at random points within the domain. Currently, it only works for problems with an analytical solution. Methodologies to obtain sensor data for problems from a file are not implemented yet. It is also not implemented for external or complex meshes.

Parameters:

Name Type Description Default
exact_solution function

The exact solution function.

 required
num_points int

The number of points to sample.

 required

Returns:

Name Type Description
Tuple

A tuple containing two arrays: sensor points and the exact solution values.
Source code in scirex/core/sciml/fe/fespace2d.py 
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
def get_sensor_data(self, exact_solution, num_points: int):
    """
    Obtain sensor data (actual solution) at random points.

    This method is used in the inverse problem to obtain the sensor data at random points within the domain.
    Currently, it only works for problems with an analytical solution.
    Methodologies to obtain sensor data for problems from a file are not implemented yet.
    It is also not implemented for external or complex meshes.

    Args:
        exact_solution (function): The exact solution function.
        num_points (int): The number of points to sample.

    Returns:
        Tuple: A tuple containing two arrays: sensor points and the exact solution values.
    """
    # generate random points within the bounds of the domain
    # get the bounds of the domain
    x_min = np.min(self.mesh.points[:, 0])
    x_max = np.max(self.mesh.points[:, 0])
    y_min = np.min(self.mesh.points[:, 1])
    y_max = np.max(self.mesh.points[:, 1])
    # sample n random points within the bounds of the domain
    # Generate points in the unit square

    num_internal_points = int(num_points * 0.9)

    points = lhs(2, samples=num_internal_points)
    points[:, 0] = x_min + (x_max - x_min) * points[:, 0]
    points[:, 1] = y_min + (y_max - y_min) * points[:, 1]
    # get the exact solution at the points
    exact_sol = exact_solution(points[:, 0], points[:, 1])

    # print the shape of the points and the exact solution
    print(f"[INFO] : Number of sensor points = {points.shape[0]}")
    print(f"[INFO] : Shape of sensor points = {points.shape}")

    # plot the points
    plt.figure(figsize=(6.4, 4.8), dpi=300)
    plt.scatter(points[:, 0], points[:, 1], marker="x", color="r", s=2)
    plt.axis("equal")
    plt.title("Sensor Points")
    plt.tight_layout()
    plt.savefig("sensor_points.png", bbox_inches="tight")

    return points, exact_sol

get_sensor_data_external(exact_sol, num_points, file_name)

This method is used to obtain the sensor data from an external file.

Parameters:

Name Type Description Default
exact_sol function

The exact solution function.

 required
num_points int

The number of points to sample.

 required
file_name str

The name of the file containing the sensor data.

 required

Returns:

Name Type Description
Tuple

A tuple containing two arrays: sensor points and the exact solution values.
Note

This method reads the sensor data from a file and samples num_points from the data. The sensor data is then returned as a tuple containing the sensor points and the exact solution values.

Source code in scirex/core/sciml/fe/fespace2d.py 
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
def get_sensor_data_external(self, exact_sol, num_points: int, file_name: str):
    """
    This method is used to obtain the sensor data from an external file.

    Args:
        exact_sol (function): The exact solution function.
        num_points (int): The number of points to sample.
        file_name (str): The name of the file containing the sensor data.

    Returns:
        Tuple: A tuple containing two arrays: sensor points and the exact solution values.

    Note:
        This method reads the sensor data from a file and samples `num_points` from the data.
        The sensor data is then returned as a tuple containing the sensor points and the exact solution values.
    """
    # use pandas to read the file
    df = pd.read_csv(file_name)

    x = df.iloc[:, 0].values
    y = df.iloc[:, 1].values
    exact_sol = df.iloc[:, 2].values

    # now sample num_points from the data
    indices = np.random.randint(0, x.shape[0], num_points)

    x = x[indices]
    y = y[indices]
    exact_sol = exact_sol[indices]

    # stack them together
    points = np.stack((x, y), axis=1)

    return points, exact_sol

get_shape_function_grad_x(cell_index)

Get the gradient of the shape function with respect to the x-coordinate.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: The actual values of the shape functions on the given cell.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace2d.py 
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
def get_shape_function_grad_x(self, cell_index: int) -> np.ndarray:
    """
    Get the gradient of the shape function with respect to the x-coordinate.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: The actual values of the shape functions on the given cell.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """
    if cell_index >= len(self.fe_cell) or cell_index < 0:
        raise ValueError(
            f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
        )

    return self.fe_cell[cell_index].basis_gradx_at_quad.copy()

get_shape_function_grad_x_ref(cell_index)

Get the gradient of the shape function with respect to the x-coordinate on the reference element.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: The actual values of the shape functions on the given cell.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace2d.py 
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
def get_shape_function_grad_x_ref(self, cell_index: int) -> np.ndarray:
    """
    Get the gradient of the shape function with respect to the x-coordinate on the reference element.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: The actual values of the shape functions on the given cell.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """
    if cell_index >= len(self.fe_cell) or cell_index < 0:
        raise ValueError(
            f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
        )

    return self.fe_cell[cell_index].basis_gradx_at_quad_ref.copy()

get_shape_function_grad_y(cell_index)

Get the gradient of the shape function with respect to y at the given cell index.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: The actual values of the shape functions on the given cell.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace2d.py 
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
def get_shape_function_grad_y(self, cell_index: int) -> np.ndarray:
    """
    Get the gradient of the shape function with respect to y at the given cell index.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: The actual values of the shape functions on the given cell.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """
    if cell_index >= len(self.fe_cell) or cell_index < 0:
        raise ValueError(
            f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
        )

    return self.fe_cell[cell_index].basis_grady_at_quad.copy()

get_shape_function_grad_y_ref(cell_index)

Get the gradient of the shape function with respect to y at the reference element.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description

np.ndarray: The actual values of the shape functions on the given cell.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Note

This function returns the gradient of the shape function with respect to y at the reference element for a given cell. The shape function gradient values are stored in the basis_grady_at_quad_ref array of the corresponding finite element cell. The cell_index parameter specifies the index of the cell for which the shape function gradient is required. If the cell_index is greater than the total number of cells, a ValueError is raised. The returned gradient values are copied from the basis_grady_at_quad_ref array to ensure immutability.

Source code in scirex/core/sciml/fe/fespace2d.py 
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
def get_shape_function_grad_y_ref(self, cell_index: int):
    """
    Get the gradient of the shape function with respect to y at the reference element.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: The actual values of the shape functions on the given cell.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.

    Note:
        This function returns the gradient of the shape function with respect to y at the reference element
        for a given cell. The shape function gradient values are stored in the `basis_grady_at_quad_ref` array
        of the corresponding finite element cell. The `cell_index` parameter specifies the index of the cell
        for which the shape function gradient is required. If the `cell_index` is greater than the total number
        of cells, a `ValueError` is raised. The returned gradient values are copied from the `basis_grady_at_quad_ref` array to ensure immutability.
    """
    if cell_index >= len(self.fe_cell) or cell_index < 0:
        raise ValueError(
            f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
        )

    return self.fe_cell[cell_index].basis_grady_at_quad_ref.copy()

get_shape_function_val(cell_index)

Get the actual values of the shape functions on a given cell.

Parameters:

Name Type Description Default
cell_index int

The index of the cell.

 required

Returns:

Type Description
ndarray

np.ndarray: The actual values of the shape functions on the given cell.

Raises:

Type Description
ValueError

If the cell_index is greater than the number of cells.
Source code in scirex/core/sciml/fe/fespace2d.py 
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
def get_shape_function_val(self, cell_index: int) -> np.ndarray:
    """
    Get the actual values of the shape functions on a given cell.

    Args:
        cell_index (int): The index of the cell.

    Returns:
        np.ndarray: The actual values of the shape functions on the given cell.

    Raises:
        ValueError: If the cell_index is greater than the number of cells.
    """
    if cell_index >= len(self.fe_cell) or cell_index < 0:
        raise ValueError(
            f"cell_index should be less than {self.n_cells} and greater than or equal to 0"
        )

    return self.fe_cell[cell_index].basis_at_quad.copy()

set_finite_elements()

Assigns the finite elements to each cell.

This method initializes the finite element objects for each cell in the mesh. It creates an instance of the FE2D_Cell class for each cell, passing the necessary parameters. The finite element objects store information about the basis functions, gradients, Jacobians, quadrature points, weights, actual coordinates, and forcing functions associated with each cell.

After initializing the finite element objects, this method prints the shape details of various matrices and updates the total number of degrees of freedom (dofs) for the entire mesh.

:return: None

Source code in scirex/core/sciml/fe/fespace2d.py 
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
def set_finite_elements(self) -> None:
    """
    Assigns the finite elements to each cell.

    This method initializes the finite element objects for each cell in the mesh.
    It creates an instance of the `FE2D_Cell` class for each cell, passing the necessary parameters.
    The finite element objects store information about the basis functions, gradients, Jacobians,
    quadrature points, weights, actual coordinates, and forcing functions associated with each cell.

    After initializing the finite element objects, this method prints the shape details of various matrices
    and updates the total number of degrees of freedom (dofs) for the entire mesh.

    :return: None
    """
    progress_bar = tqdm(
        total=self.n_cells,
        desc="Fe2D_cell Setup",
        unit="cells_assembled",
        bar_format="{l_bar}{bar:40}{r_bar}{bar:-10b}",
        colour="blue",
        ncols=100,
    )

    dof = 0
    for i in range(self.n_cells):
        self.fe_cell.append(
            FE2D_Cell(
                self.cells[i],
                self.cell_type,
                self.fe_order,
                self.fe_type,
                self.quad_order,
                self.quad_type,
                self.fe_transformation_type,
                self.forcing_function,
            )
        )

        # obtain the shape of the basis function (n_test, N_quad)
        dof += self.fe_cell[i].basis_at_quad.shape[1]

        progress_bar.update(1)
    # print the Shape details of all the matrices from cell 0 using print_table function
    title = [
        "Shape function Matrix Shape",
        "Shape function Gradient Matrix Shape",
        "Jacobian Matrix Shape",
        "Quadrature Points Shape",
        "Quadrature Weights Shape",
        "Quadrature Actual Coordinates Shape",
        "Forcing Function Shape",
    ]
    values = [
        self.fe_cell[0].basis_at_quad.shape,
        self.fe_cell[0].basis_gradx_at_quad.shape,
        self.fe_cell[0].jacobian.shape,
        self.fe_cell[0].quad_xi.shape,
        self.fe_cell[0].quad_weight.shape,
        self.fe_cell[0].quad_actual_coordinates.shape,
        self.fe_cell[0].forcing_at_quad.shape,
    ]
    print_table("fe Matrix Shapes", ["Matrix", "Shape"], title, values)

    # update the total number of dofs
    self.total_dofs = dof

Geometry

Abstract base class for geometry and mesh operations.

This class defines the interface that all geometry implementations must follow, providing the basic structure for mesh handling operations in both 2D and 3D contexts.

Attributes:

Name Type Description
mesh_type

Type of mesh elements (e.g., 'quadrilateral', 'triangle')
mesh_generation_method

Method for mesh generation ('internal'/'external')
Example

class Geometry2D(Geometry): ... def init(self, mesh_type='quadrilateral', ... method='internal'): ... super().init(mesh_type, method) ... ... def read_mesh(self, mesh_file, boundary_level, ... sampling_method, refine_level): ... # Implementation ... pass ... ... def generate_vtk_for_test(self): ... # Implementation ... pass ... ... def get_test_points(self): ... # Implementation ... return points

Note

This is an abstract base class. Concrete implementations must override: - read_mesh() - generate_vtk_for_test() - get_test_points()

Each implementation should provide appropriate mesh handling for its specific dimensional and element type requirements.

Source code in scirex/core/sciml/geometry/geometry.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
class Geometry:
    """Abstract base class for geometry and mesh operations.

    This class defines the interface that all geometry implementations must
    follow, providing the basic structure for mesh handling operations in
    both 2D and 3D contexts.

    Attributes:
        mesh_type: Type of mesh elements (e.g., 'quadrilateral', 'triangle')
        mesh_generation_method: Method for mesh generation ('internal'/'external')

    Example:
        >>> class Geometry2D(Geometry):
        ...     def __init__(self, mesh_type='quadrilateral',
        ...                  method='internal'):
        ...         super().__init__(mesh_type, method)
        ...
        ...     def read_mesh(self, mesh_file, boundary_level,
        ...                   sampling_method, refine_level):
        ...         # Implementation
        ...         pass
        ...
        ...     def generate_vtk_for_test(self):
        ...         # Implementation
        ...         pass
        ...
        ...     def get_test_points(self):
        ...         # Implementation
        ...         return points

    Note:
        This is an abstract base class. Concrete implementations must override:
        - read_mesh()
        - generate_vtk_for_test()
        - get_test_points()

        Each implementation should provide appropriate mesh handling for its
        specific dimensional and element type requirements.
    """

    def __init__(self, mesh_type: str, mesh_generation_method: str):
        """
        Constructor for the Geometry class.

        Args:
            mesh_type: Type of mesh elements (e.g., 'quadrilateral', 'triangle')
            mesh_generation_method: Method for mesh generation ('internal'/'external')

        Returns:
            None
        """
        self.mesh_type = mesh_type
        self.mesh_generation_method = mesh_generation_method

    @abstractmethod
    def read_mesh(
        self,
        mesh_file: str,
        boundary_point_refinement_level: int,
        bd_sampling_method: str,
        refinement_level: int,
    ):
        """
        Abstract method to read mesh from Gmsh. This method should be implemented by the derived classes.

        Args:
            mesh_file (str): Path to the mesh file
            boundary_point_refinement_level (int): Level of refinement for boundary points
            bd_sampling_method (str): Sampling method for boundary points
            refinement_level (int): Level of mesh refinement

        Returns:
            None
        """

    @abstractmethod
    def generate_vtk_for_test(self):
        """
        Generates a VTK from Mesh file (External) or using gmsh (for Internal).

        Args:
        None

        Returns:
        None
        """

    @abstractmethod
    def get_test_points(self):
        """
        This function is used to extract the test points from the given mesh

        Args:
            None

        Returns:
            points (np.ndarray): Test points extracted from the mesh
        """

__init__(mesh_type, mesh_generation_method)

Constructor for the Geometry class.

Parameters:

Name Type Description Default
mesh_type str

Type of mesh elements (e.g., 'quadrilateral', 'triangle')

 required
mesh_generation_method str

Method for mesh generation ('internal'/'external')

 required

Returns:

Type Description

None
Source code in scirex/core/sciml/geometry/geometry.py 
113
114
115
116
117
118
119
120
121
122
123
124
125
def __init__(self, mesh_type: str, mesh_generation_method: str):
    """
    Constructor for the Geometry class.

    Args:
        mesh_type: Type of mesh elements (e.g., 'quadrilateral', 'triangle')
        mesh_generation_method: Method for mesh generation ('internal'/'external')

    Returns:
        None
    """
    self.mesh_type = mesh_type
    self.mesh_generation_method = mesh_generation_method

generate_vtk_for_test() abstractmethod

Generates a VTK from Mesh file (External) or using gmsh (for Internal).

Args: None

Returns: None

Source code in scirex/core/sciml/geometry/geometry.py 
148
149
150
151
152
153
154
155
156
157
158
@abstractmethod
def generate_vtk_for_test(self):
    """
    Generates a VTK from Mesh file (External) or using gmsh (for Internal).

    Args:
    None

    Returns:
    None
    """

get_test_points() abstractmethod

This function is used to extract the test points from the given mesh

Returns:

Name Type Description
points ndarray

Test points extracted from the mesh
Source code in scirex/core/sciml/geometry/geometry.py 
160
161
162
163
164
165
166
167
168
169
170
@abstractmethod
def get_test_points(self):
    """
    This function is used to extract the test points from the given mesh

    Args:
        None

    Returns:
        points (np.ndarray): Test points extracted from the mesh
    """

read_mesh(mesh_file, boundary_point_refinement_level, bd_sampling_method, refinement_level) abstractmethod

Abstract method to read mesh from Gmsh. This method should be implemented by the derived classes.

Parameters:

Name Type Description Default
mesh_file str

Path to the mesh file

 required
boundary_point_refinement_level int

Level of refinement for boundary points

 required
bd_sampling_method str

Sampling method for boundary points

 required
refinement_level int

Level of mesh refinement

 required

Returns:

Type Description

None
Source code in scirex/core/sciml/geometry/geometry.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
@abstractmethod
def read_mesh(
    self,
    mesh_file: str,
    boundary_point_refinement_level: int,
    bd_sampling_method: str,
    refinement_level: int,
):
    """
    Abstract method to read mesh from Gmsh. This method should be implemented by the derived classes.

    Args:
        mesh_file (str): Path to the mesh file
        boundary_point_refinement_level (int): Level of refinement for boundary points
        bd_sampling_method (str): Sampling method for boundary points
        refinement_level (int): Level of mesh refinement

    Returns:
        None
    """

Geometry_2D

Bases: Geometry

Implements 2D geometry and mesh handling capabilities.

This class provides comprehensive functionality for managing 2D meshes, including both internal generation and external mesh reading. It supports various mesh operations, boundary handling, and visualization capabilities.

Attributes:

Name Type Description
mesh_type

Type of mesh elements ('quadrilateral')
mesh_generation_method

Method of mesh generation ('internal'/'external')
n_test_points_x

Number of test points in x-direction
n_test_points_y

Number of test points in y-direction
output_folder

Path for output files
is_optimized

Flag for mesh optimization
n_cells_x

Number of cells in x-direction (internal mesh)
n_cells_y

Number of cells in y-direction (internal mesh)
x_limits

Domain limits in x-direction
y_limits

Domain limits in y-direction
mesh_file_name

Name of external mesh file
mesh

MeshIO mesh object
bd_dict

Dictionary of boundary points
cell_points

Array of cell vertices
test_points

Array of test points
Example

geometry = Geometry_2D( ... mesh_type='quadrilateral', ... mesh_generation_method='internal', ... n_test_points_x=10, ... n_test_points_y=10, ... output_folder='./output' ... ) cells, bounds = geometry.generate_quad_mesh_internal( ... x_limits=(0,1), ... y_limits=(0,1), ... n_cells_x=5, ... n_cells_y=5, ... num_boundary_points=40 ... )

Note
  • Only supports quadrilateral elements
  • Internal mesh generation is limited to rectangular domains
  • External mesh reading requires Gmsh format
  • Boundary points can be sampled uniformly or using LHS
Source code in scirex/core/sciml/geometry/geometry_2d.py 
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
class Geometry_2D(Geometry):
    """Implements 2D geometry and mesh handling capabilities.

    This class provides comprehensive functionality for managing 2D meshes,
    including both internal generation and external mesh reading. It supports
    various mesh operations, boundary handling, and visualization capabilities.

    Attributes:
        mesh_type: Type of mesh elements ('quadrilateral')
        mesh_generation_method: Method of mesh generation ('internal'/'external')
        n_test_points_x: Number of test points in x-direction
        n_test_points_y: Number of test points in y-direction
        output_folder: Path for output files
        is_optimized: Flag for mesh optimization
        n_cells_x: Number of cells in x-direction (internal mesh)
        n_cells_y: Number of cells in y-direction (internal mesh)
        x_limits: Domain limits in x-direction
        y_limits: Domain limits in y-direction
        mesh_file_name: Name of external mesh file
        mesh: MeshIO mesh object
        bd_dict: Dictionary of boundary points
        cell_points: Array of cell vertices
        test_points: Array of test points

    Example:
        >>> geometry = Geometry_2D(
        ...     mesh_type='quadrilateral',
        ...     mesh_generation_method='internal',
        ...     n_test_points_x=10,
        ...     n_test_points_y=10,
        ...     output_folder='./output'
        ... )
        >>> cells, bounds = geometry.generate_quad_mesh_internal(
        ...     x_limits=(0,1),
        ...     y_limits=(0,1),
        ...     n_cells_x=5,
        ...     n_cells_y=5,
        ...     num_boundary_points=40
        ... )

    Note:
        - Only supports quadrilateral elements
        - Internal mesh generation is limited to rectangular domains
        - External mesh reading requires Gmsh format
        - Boundary points can be sampled uniformly or using LHS
    """

    def __init__(
        self,
        mesh_type: str,
        mesh_generation_method: str,
        n_test_points_x: int,
        n_test_points_y: int,
        output_folder: str,
        is_optimized: bool = False,
    ):
        """
        Constructor for Geometry_2D class.

        Args:
            mesh_type: Type of mesh elements ('quadrilateral')
            mesh_generation_method: Method of mesh generation ('internal'/'external')
            n_test_points_x: Number of test points in x-direction
            n_test_points_y: Number of test points in y-direction
            output_folder: Path for output files
            is_optimized: Flag for mesh optimization

        Raises:
            ValueError: If mesh type or generation method is invalid

        Returns:
            None
        """
        # Call the super class constructor
        super().__init__(mesh_type, mesh_generation_method)
        self.mesh_type = mesh_type
        self.mesh_generation_method = mesh_generation_method
        self.n_test_points_x = n_test_points_x
        self.n_test_points_y = n_test_points_y
        self.output_folder = output_folder
        self.is_optimized = is_optimized

        if self.mesh_generation_method not in ["internal", "external"]:
            print(
                f"Invalid mesh generation method {self.mesh_generation_method} in {self.__class__.__name__} from {__name__}."
            )
            raise ValueError(
                "Mesh generation method should be either internal or external."
            )

        if self.mesh_type not in ["quadrilateral"]:
            print(
                f"Invalid mesh type {self.mesh_type} in {self.__class__.__name__} from {__name__}."
            )
            raise ValueError("Mesh type should be quadrilateral only.")

        # To be filled - only when mesh is internal
        self.n_cells_x = None
        self.n_cells_y = None
        self.x_limits = None
        self.y_limits = None

        # to be filled by external
        self.mesh_file_name = None
        self.mesh = None
        self.bd_dict = None
        self.cell_points = None
        self.test_points = None

    def read_mesh(
        self,
        mesh_file: str,
        boundary_point_refinement_level: int,
        bd_sampling_method: str,
        refinement_level: int,
    ):
        """
        Reads mesh from a Gmsh .msh file and extracts cell information.

        Args:
            mesh_file: Path to the mesh file
            boundary_point_refinement_level: Level of boundary point refinement
            bd_sampling_method: Method for boundary point sampling ('uniform'/'lhs')
            refinement_level: Level of mesh refinement

        Returns:
            cell_points: Array of cell vertices
            bd_dict: Dictionary of boundary points

        Raises:
            ValueError: If mesh file format is invalid
        """

        self.mesh_file_name = mesh_file

        # bd_sampling_method = "uniform"  # "uniform" or "lhs"

        file_extension = Path(mesh_file).suffix

        if file_extension != ".mesh":
            raise ValueError("Mesh file should be in .mesh format.")

        # Read mesh using meshio
        self.mesh = meshio.read(mesh_file)

        if self.mesh_type == "quadrilateral":
            # Extract cell information
            cells = self.mesh.cells_dict["quad"]

        num_cells = cells.shape[0]
        print(f"[INFO] : Number of cells = {num_cells}")
        cell_points = self.mesh.points[cells][
            :, :, 0:2
        ]  # remove the z coordinate, which is 0 for all points

        # loop over all cells and rearrange the points in anticlockwise direction
        for i in range(num_cells):
            cell = cell_points[i]
            # get the centroid of the cell
            centroid = np.mean(cell, axis=0)
            # get the angle of each point with respect to the centroid
            angles = np.arctan2(cell[:, 1] - centroid[1], cell[:, 0] - centroid[0])
            # sort the points based on the angles
            cell_points[i] = cell[np.argsort(angles)]

        # Extract number of points within each cell
        print(f"[INFO] : Number of points per cell = {cell_points.shape}")

        # Collect the Boundary point id's within the domain
        boundary_edges = self.mesh.cells_dict["line"]

        # Using the point id, collect the coordinates of the boundary points
        boundary_coordinates = self.mesh.points[boundary_edges]

        # Number of Existing Boundary points
        print(
            f"[INFO] : Number of Bound points before refinement = {np.unique(boundary_coordinates.reshape(-1,3)).shape[0] * 0.5 + 1}"
        )

        # now Get the physical tag of the boundary edges
        boundary_tags = self.mesh.cell_data["medit:ref"][0]

        # Generate a Dictionary of boundary tags and boundary coordinates
        # Keys will be the boundary tags and values will be the list of coordinates
        boundary_dict = {}

        # refine the boundary points based on the number of boundary points needed
        for i in range(boundary_coordinates.shape[0]):
            p1 = boundary_coordinates[i, 0, :]
            p2 = boundary_coordinates[i, 1, :]

            if bd_sampling_method == "uniform":
                # take the current point and next point and then perform a uniform sampling
                new_points = np.linspace(
                    p1, p2, pow(2, boundary_point_refinement_level) + 1
                )
            elif bd_sampling_method == "lhs":
                # take the current point and next point and then perform a uniform sampling
                new_points = lhs(2, pow(2, boundary_point_refinement_level) + 1)
                new_points[:, 0] = new_points[:, 0] * (p2[0] - p1[0]) + p1[0]
                new_points[:, 1] = new_points[:, 1] * (p2[1] - p1[1]) + p1[1]
            else:
                print(
                    f"Invalid sampling method {bd_sampling_method} in {self.__class__.__name__} from {__name__}."
                )
                raise ValueError("Sampling method should be either uniform or lhs.")

            # get the boundary tag
            tag = boundary_tags[i]

            if tag not in boundary_dict:
                boundary_dict[tag] = new_points
            else:
                current_val = new_points
                prev_val = boundary_dict[tag]
                final = np.vstack([prev_val, current_val])
                boundary_dict[tag] = final

        # get unique
        for tag in boundary_dict.keys():
            val = boundary_dict[tag]
            val = np.unique(val, axis=0)
            boundary_dict[tag] = val

        self.bd_dict = boundary_dict
        # print the new boundary points  on each boundary tag (key) in a tabular format

        total_bound_points = 0
        print(f"| {'Boundary ID':<12} | {'Number of Points':<16} |")
        print(f"| {'-'*12:<12}---{'-'*16:<16} |")
        for k, v in self.bd_dict.items():
            print(f"| {k:<12} | {v.shape[0]:<16} |")
            total_bound_points += v.shape[0]

        print(f"[INFO] : No of bound pts after refinement:  {total_bound_points}")

        # Assign to class values
        self.cell_points = cell_points

        # generate testvtk
        self.generate_vtk_for_test()

        return cell_points, self.bd_dict

    def generate_quad_mesh_internal(
        self,
        x_limits: tuple,
        y_limits: tuple,
        n_cells_x: int,
        n_cells_y: int,
        num_boundary_points: int,
    ):
        """
        Generate and save a quadrilateral mesh with physical curves.

        Args:
            x_limits: Domain limits in x-direction
            y_limits: Domain limits in y-direction
            n_cells_x: Number of cells in x-direction
            n_cells_y: Number of cells in y-direction
            num_boundary_points: Number of boundary points

        Returns:
            cell_points: Array of cell vertices
            bd_dict: Dictionary of boundary points
        """

        self.n_cells_x = n_cells_x
        self.n_cells_y = n_cells_y
        self.x_limits = x_limits
        self.y_limits = y_limits

        # generate linspace of points in x and y direction
        x = np.linspace(x_limits[0], x_limits[1], n_cells_x + 1)
        y = np.linspace(y_limits[0], y_limits[1], n_cells_y + 1)

        # Generate quad cells from the points
        # the output should be a list of 4 points for each cell , each being a list of 2 points [x,y]
        cells = []

        for i in range(n_cells_x):
            for j in range(n_cells_y):
                # get the four points of the cell
                p1 = [x[i], y[j]]
                p2 = [x[i + 1], y[j]]
                p3 = [x[i + 1], y[j + 1]]
                p4 = [x[i], y[j + 1]]

                # append the points to the cells
                cells.append([p1, p2, p3, p4])

        # convert to numpy array
        cells = np.array(cells, dtype=np.float64)

        # use arctan2 to sort the points in anticlockwise direction
        # loop over all cells and rearrange the points in anticlockwise direction
        for i in range(cells.shape[0]):
            cell = cells[i]
            # get the centroid of the cell
            centroid = np.mean(cell, axis=0)
            # get the angle of each point with respect to the centroid
            angles = np.arctan2(cell[:, 1] - centroid[1], cell[:, 0] - centroid[0])
            # sort the points based on the angles
            cells[i] = cell[np.argsort(angles)]

        # generate a meshio mesh object using the cells
        self.mesh = meshio.Mesh(
            points=cells.reshape(-1, 2), cells=[("quad", cells.reshape(-1, 4))]
        )

        # lets generate the boundary points, this function will return a dictionary of boundary points
        # the keys will be the boundary tags and values will be the list of boundary points
        bd_points = {}

        num_bound_per_side = int(num_boundary_points / 4)

        def _temp_bd_func(start, end, num_pts):
            """
            This function returns the boundary points between the start and end points
            using lhs sampling.

            Args:
                start: Start point of the boundary
                end: End point of the boundary
                num_pts: Number of boundary points to be generated

            Returns:
                bd_pts: Array of boundary points
            """
            # generate the boundary points using lhs as a np.float64 array
            bd_pts = lhs(1, num_pts).astype(np.float64)
            # scale the points
            bd_pts = bd_pts * (end - start) + start

            return bd_pts.reshape(-1)

        # bottom boundary
        y_bottom = (
            np.ones(num_bound_per_side, dtype=np.float64) * y_limits[0]
        ).reshape(-1)
        x_bottom = _temp_bd_func(x_limits[0], x_limits[1], num_bound_per_side)
        bd_points[1000] = np.vstack([x_bottom, y_bottom]).T

        # right boundary
        x_right = (np.ones(num_bound_per_side, dtype=np.float64) * x_limits[1]).reshape(
            -1
        )
        y_right = _temp_bd_func(y_limits[0], y_limits[1], num_bound_per_side)
        bd_points[1001] = np.vstack([x_right, y_right]).T

        # top boundary
        y_top = (np.ones(num_bound_per_side, dtype=np.float64) * y_limits[1]).reshape(
            -1
        )
        x_top = _temp_bd_func(x_limits[0], x_limits[1], num_bound_per_side)
        bd_points[1002] = np.vstack([x_top, y_top]).T

        # left boundary
        x_left = (np.ones(num_bound_per_side, dtype=np.float64) * x_limits[0]).reshape(
            -1
        )
        y_left = _temp_bd_func(y_limits[0], y_limits[1], num_bound_per_side)
        bd_points[1003] = np.vstack([x_left, y_left]).T

        self.cell_points = cells
        self.bd_dict = bd_points

        # generate vtk
        if not self.is_optimized:
            self.generate_vtk_for_test()

        return self.cell_points, self.bd_dict

    def generate_vtk_for_test(self):
        """
        Generates a VTK from Mesh file (External) or using gmsh (for Internal).

        Args:
            None

        Returns:
            None
        """

        if self.mesh_generation_method == "internal":
            # initialise the mesh
            gmsh.initialize()

            # Now, lets generate the mesh with the points.
            x_range = self.x_limits[1] - self.x_limits[0]
            y_range = self.y_limits[1] - self.y_limits[0]

            mesh_size_x = x_range / self.n_test_points_x
            mesh_size_y = y_range / self.n_test_points_y

            # generate a gmsh with the given parameters
            Xmin = self.x_limits[0]
            Xmax = self.x_limits[1]
            Ymin = self.y_limits[0]
            Ymax = self.y_limits[1]

            point1 = gmsh.model.geo.add_point(Xmin, Ymin, 0, mesh_size_x)
            point2 = gmsh.model.geo.add_point(Xmax, Ymin, 0, mesh_size_x)
            point3 = gmsh.model.geo.add_point(Xmax, Ymax, 0, mesh_size_y)
            point4 = gmsh.model.geo.add_point(Xmin, Ymax, 0, mesh_size_y)

            line1 = gmsh.model.geo.add_line(point1, point2, 1000)  ## Bottom
            line2 = gmsh.model.geo.add_line(point2, point3, 1001)  ## Right
            line3 = gmsh.model.geo.add_line(point3, point4, 1002)  ## Top
            line4 = gmsh.model.geo.add_line(point4, point1, 1003)  ## Left

            face1 = gmsh.model.geo.add_curve_loop([line1, line2, line3, line4])

            gmsh.model.geo.add_plane_surface([face1])

            # Create the relevant Gmsh data structures
            # from Gmsh model.
            gmsh.model.geo.synchronize()

            # Generate mesh:
            gmsh.model.mesh.generate()

            mesh_file_name = Path(self.output_folder) / "internal.msh"
            vtk_file_name = Path(self.output_folder) / "internal.vtk"

            gmsh.write(str(mesh_file_name))
            print("[INFO] : Internal mesh file generated at ", str(mesh_file_name))

            # close the gmsh
            gmsh.finalize()

            # read the mesh using meshio
            mesh = meshio.gmsh.read(str(mesh_file_name))
            meshio.vtk.write(str(vtk_file_name), mesh, binary=False, fmt_version="4.2")

            print(
                "[INFO] : VTK file for internal mesh file generated at ",
                str(mesh_file_name),
            )

        elif self.mesh_generation_method == "external":

            vtk_file_name = Path(self.output_folder) / "external.vtk"

            # Use the internal mesh to generate the vtk file
            mesh = meshio.read(str(self.mesh_file_name))
            meshio.vtk.write(str(vtk_file_name), mesh, binary=False, fmt_version="4.2")

            print(
                "[INFO] : VTK file for external mesh file generated at ",
                str(vtk_file_name),
            )

    def get_test_points(self):
        """
        This function is used to extract the test points from the given mesh

        Args:
            None

        Returns:
            test_points (np.ndarray): Array of test points
        """

        if self.mesh_generation_method == "internal":
            # vtk_file_name  = Path(self.output_folder) / "internal.vtk"
            # code over written to plot from np.linspace instead of vtk file
            # generate linspace of points in x and y direction based on x and y limits
            x = np.linspace(self.x_limits[0], self.x_limits[1], self.n_test_points_x)
            y = np.linspace(self.y_limits[0], self.y_limits[1], self.n_test_points_y)
            # generate meshgrid
            x_grid, y_grid = np.meshgrid(x, y)
            # stack the points
            self.test_points = np.vstack([x_grid.flatten(), y_grid.flatten()]).T

            return self.test_points

        elif self.mesh_generation_method == "external":
            vtk_file_name = Path(self.output_folder) / "external.vtk"

        mesh = meshio.read(str(vtk_file_name))
        points = mesh.points
        return points[:, 0:2]  # return only first two columns

    def write_vtk(
        self, solution: np.ndarray, output_path: str, filename: str, data_names: list
    ):
        """
        Writes the data to a VTK file.

        Args:
            solution: The solution data to be written
            output_path: The output path for the VTK file
            filename: The name of the output file
            data_names: List of data names

        Returns:
            None
        """
        # read the existing vtk into file
        if self.mesh_generation_method == "internal":
            vtk_file_name = Path(self.output_folder) / "internal.vtk"
        elif self.mesh_generation_method == "external":
            vtk_file_name = Path(self.output_folder) / "external.vtk"

        data = []
        with open(vtk_file_name, "r", encoding="utf-8") as File:
            for line in File:
                data.append(line)

        # get the output file name
        output_file_name = Path(output_path) / filename

        if solution.shape[1] != len(data_names):
            print("[Error] : File : geometry_2d.py, Function: write_vtk")
            print(
                "Num Columns in solution = ",
                solution.shape[1],
                " Num of data names = ",
                len(data_names),
            )
            raise ValueError("Number of data names and solution columns are not equal")

        # write the data to the output file
        with open(str(output_file_name), "w", encoding="utf-8") as FN:
            for line in data:
                FN.write(line)
                if "POINT_DATA" in line.strip():
                    break

            for i in range(solution.shape[1]):
                FN.write("SCALARS " + data_names[i] + " float\n")
                FN.write("LOOKUP_TABLE default\n")
                np.savetxt(FN, solution[:, i])
                FN.write("\n")

        # save the vtk file as image
        # self.save_vtk_as_image(str(output_file_name), data_names)

    def plot_adaptive_mesh(
        self, cells_list, area_averaged_cell_loss_list, epoch, filename="cell_residual"
    ):
        """
        Plots the residuals in each cell of the mesh.

        Args:
            cells_list: List of cell vertices
            area_averaged_cell_loss_list: List of area averaged cell loss
            epoch: The epoch number
            filename: The output filename

        Returns:
            None
        """

        plt.figure(figsize=(6.4, 4.8), dpi=300)

        # normalise colors
        norm = mcolors.Normalize(
            vmin=np.min(area_averaged_cell_loss_list),
            vmax=np.max(area_averaged_cell_loss_list),
        )

        # Create a colormap
        colormap = plt.cm.jet

        for index, cell in enumerate(cells_list):
            x = cell[:, 0]
            y = cell[:, 1]

            x = np.append(x, x[0])
            y = np.append(y, y[0])

            curr_cell_loss = float(area_averaged_cell_loss_list[index])

            color = colormap(norm(curr_cell_loss))

            plt.fill(x, y, color=color, alpha=0.9)

            plt.plot(x, y, "k")

            # # compute x_min, x_max, y_min, y_max
            # x_min = np.min(x)
            # x_max = np.max(x)
            # y_min = np.min(y)
            # y_max = np.max(y)

            # # compute centroid of the cells
            # centroid = np.array([np.mean(x), np.mean(y)])

            # plot the loss text within the cell
            # plt.text(centroid[0], centroid[1], f"{curr_cell_loss:.3e}", fontsize=16, horizontalalignment='center', verticalalignment='center')

        sm = plt.cm.ScalarMappable(cmap=colormap, norm=norm)
        sm.set_array([])
        plt.colorbar(sm)

        # output filename
        output_filename = Path(f"{self.output_folder}/{filename}_{epoch}.png")
        plt.title(f"Cell Residual")
        plt.savefig(str(output_filename), dpi=300)

__init__(mesh_type, mesh_generation_method, n_test_points_x, n_test_points_y, output_folder, is_optimized=False)

Constructor for Geometry_2D class.

Parameters:

Name Type Description Default
mesh_type str

Type of mesh elements ('quadrilateral')

 required
mesh_generation_method str

Method of mesh generation ('internal'/'external')

 required
n_test_points_x int

Number of test points in x-direction

 required
n_test_points_y int

Number of test points in y-direction

 required
output_folder str

Path for output files

 required
is_optimized bool

Flag for mesh optimization

 False

Raises:

Type Description
ValueError

If mesh type or generation method is invalid

Returns:

Type Description

None
Source code in scirex/core/sciml/geometry/geometry_2d.py 
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
def __init__(
    self,
    mesh_type: str,
    mesh_generation_method: str,
    n_test_points_x: int,
    n_test_points_y: int,
    output_folder: str,
    is_optimized: bool = False,
):
    """
    Constructor for Geometry_2D class.

    Args:
        mesh_type: Type of mesh elements ('quadrilateral')
        mesh_generation_method: Method of mesh generation ('internal'/'external')
        n_test_points_x: Number of test points in x-direction
        n_test_points_y: Number of test points in y-direction
        output_folder: Path for output files
        is_optimized: Flag for mesh optimization

    Raises:
        ValueError: If mesh type or generation method is invalid

    Returns:
        None
    """
    # Call the super class constructor
    super().__init__(mesh_type, mesh_generation_method)
    self.mesh_type = mesh_type
    self.mesh_generation_method = mesh_generation_method
    self.n_test_points_x = n_test_points_x
    self.n_test_points_y = n_test_points_y
    self.output_folder = output_folder
    self.is_optimized = is_optimized

    if self.mesh_generation_method not in ["internal", "external"]:
        print(
            f"Invalid mesh generation method {self.mesh_generation_method} in {self.__class__.__name__} from {__name__}."
        )
        raise ValueError(
            "Mesh generation method should be either internal or external."
        )

    if self.mesh_type not in ["quadrilateral"]:
        print(
            f"Invalid mesh type {self.mesh_type} in {self.__class__.__name__} from {__name__}."
        )
        raise ValueError("Mesh type should be quadrilateral only.")

    # To be filled - only when mesh is internal
    self.n_cells_x = None
    self.n_cells_y = None
    self.x_limits = None
    self.y_limits = None

    # to be filled by external
    self.mesh_file_name = None
    self.mesh = None
    self.bd_dict = None
    self.cell_points = None
    self.test_points = None

generate_quad_mesh_internal(x_limits, y_limits, n_cells_x, n_cells_y, num_boundary_points)

Generate and save a quadrilateral mesh with physical curves.

Parameters:

Name Type Description Default
x_limits tuple

Domain limits in x-direction

 required
y_limits tuple

Domain limits in y-direction

 required
n_cells_x int

Number of cells in x-direction

 required
n_cells_y int

Number of cells in y-direction

 required
num_boundary_points int

Number of boundary points

 required

Returns:

Name Type Description
cell_points

Array of cell vertices
bd_dict

Dictionary of boundary points
Source code in scirex/core/sciml/geometry/geometry_2d.py 
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
def generate_quad_mesh_internal(
    self,
    x_limits: tuple,
    y_limits: tuple,
    n_cells_x: int,
    n_cells_y: int,
    num_boundary_points: int,
):
    """
    Generate and save a quadrilateral mesh with physical curves.

    Args:
        x_limits: Domain limits in x-direction
        y_limits: Domain limits in y-direction
        n_cells_x: Number of cells in x-direction
        n_cells_y: Number of cells in y-direction
        num_boundary_points: Number of boundary points

    Returns:
        cell_points: Array of cell vertices
        bd_dict: Dictionary of boundary points
    """

    self.n_cells_x = n_cells_x
    self.n_cells_y = n_cells_y
    self.x_limits = x_limits
    self.y_limits = y_limits

    # generate linspace of points in x and y direction
    x = np.linspace(x_limits[0], x_limits[1], n_cells_x + 1)
    y = np.linspace(y_limits[0], y_limits[1], n_cells_y + 1)

    # Generate quad cells from the points
    # the output should be a list of 4 points for each cell , each being a list of 2 points [x,y]
    cells = []

    for i in range(n_cells_x):
        for j in range(n_cells_y):
            # get the four points of the cell
            p1 = [x[i], y[j]]
            p2 = [x[i + 1], y[j]]
            p3 = [x[i + 1], y[j + 1]]
            p4 = [x[i], y[j + 1]]

            # append the points to the cells
            cells.append([p1, p2, p3, p4])

    # convert to numpy array
    cells = np.array(cells, dtype=np.float64)

    # use arctan2 to sort the points in anticlockwise direction
    # loop over all cells and rearrange the points in anticlockwise direction
    for i in range(cells.shape[0]):
        cell = cells[i]
        # get the centroid of the cell
        centroid = np.mean(cell, axis=0)
        # get the angle of each point with respect to the centroid
        angles = np.arctan2(cell[:, 1] - centroid[1], cell[:, 0] - centroid[0])
        # sort the points based on the angles
        cells[i] = cell[np.argsort(angles)]

    # generate a meshio mesh object using the cells
    self.mesh = meshio.Mesh(
        points=cells.reshape(-1, 2), cells=[("quad", cells.reshape(-1, 4))]
    )

    # lets generate the boundary points, this function will return a dictionary of boundary points
    # the keys will be the boundary tags and values will be the list of boundary points
    bd_points = {}

    num_bound_per_side = int(num_boundary_points / 4)

    def _temp_bd_func(start, end, num_pts):
        """
        This function returns the boundary points between the start and end points
        using lhs sampling.

        Args:
            start: Start point of the boundary
            end: End point of the boundary
            num_pts: Number of boundary points to be generated

        Returns:
            bd_pts: Array of boundary points
        """
        # generate the boundary points using lhs as a np.float64 array
        bd_pts = lhs(1, num_pts).astype(np.float64)
        # scale the points
        bd_pts = bd_pts * (end - start) + start

        return bd_pts.reshape(-1)

    # bottom boundary
    y_bottom = (
        np.ones(num_bound_per_side, dtype=np.float64) * y_limits[0]
    ).reshape(-1)
    x_bottom = _temp_bd_func(x_limits[0], x_limits[1], num_bound_per_side)
    bd_points[1000] = np.vstack([x_bottom, y_bottom]).T

    # right boundary
    x_right = (np.ones(num_bound_per_side, dtype=np.float64) * x_limits[1]).reshape(
        -1
    )
    y_right = _temp_bd_func(y_limits[0], y_limits[1], num_bound_per_side)
    bd_points[1001] = np.vstack([x_right, y_right]).T

    # top boundary
    y_top = (np.ones(num_bound_per_side, dtype=np.float64) * y_limits[1]).reshape(
        -1
    )
    x_top = _temp_bd_func(x_limits[0], x_limits[1], num_bound_per_side)
    bd_points[1002] = np.vstack([x_top, y_top]).T

    # left boundary
    x_left = (np.ones(num_bound_per_side, dtype=np.float64) * x_limits[0]).reshape(
        -1
    )
    y_left = _temp_bd_func(y_limits[0], y_limits[1], num_bound_per_side)
    bd_points[1003] = np.vstack([x_left, y_left]).T

    self.cell_points = cells
    self.bd_dict = bd_points

    # generate vtk
    if not self.is_optimized:
        self.generate_vtk_for_test()

    return self.cell_points, self.bd_dict

generate_vtk_for_test()

Generates a VTK from Mesh file (External) or using gmsh (for Internal).

Returns:

Type Description

None
Source code in scirex/core/sciml/geometry/geometry_2d.py 
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
def generate_vtk_for_test(self):
    """
    Generates a VTK from Mesh file (External) or using gmsh (for Internal).

    Args:
        None

    Returns:
        None
    """

    if self.mesh_generation_method == "internal":
        # initialise the mesh
        gmsh.initialize()

        # Now, lets generate the mesh with the points.
        x_range = self.x_limits[1] - self.x_limits[0]
        y_range = self.y_limits[1] - self.y_limits[0]

        mesh_size_x = x_range / self.n_test_points_x
        mesh_size_y = y_range / self.n_test_points_y

        # generate a gmsh with the given parameters
        Xmin = self.x_limits[0]
        Xmax = self.x_limits[1]
        Ymin = self.y_limits[0]
        Ymax = self.y_limits[1]

        point1 = gmsh.model.geo.add_point(Xmin, Ymin, 0, mesh_size_x)
        point2 = gmsh.model.geo.add_point(Xmax, Ymin, 0, mesh_size_x)
        point3 = gmsh.model.geo.add_point(Xmax, Ymax, 0, mesh_size_y)
        point4 = gmsh.model.geo.add_point(Xmin, Ymax, 0, mesh_size_y)

        line1 = gmsh.model.geo.add_line(point1, point2, 1000)  ## Bottom
        line2 = gmsh.model.geo.add_line(point2, point3, 1001)  ## Right
        line3 = gmsh.model.geo.add_line(point3, point4, 1002)  ## Top
        line4 = gmsh.model.geo.add_line(point4, point1, 1003)  ## Left

        face1 = gmsh.model.geo.add_curve_loop([line1, line2, line3, line4])

        gmsh.model.geo.add_plane_surface([face1])

        # Create the relevant Gmsh data structures
        # from Gmsh model.
        gmsh.model.geo.synchronize()

        # Generate mesh:
        gmsh.model.mesh.generate()

        mesh_file_name = Path(self.output_folder) / "internal.msh"
        vtk_file_name = Path(self.output_folder) / "internal.vtk"

        gmsh.write(str(mesh_file_name))
        print("[INFO] : Internal mesh file generated at ", str(mesh_file_name))

        # close the gmsh
        gmsh.finalize()

        # read the mesh using meshio
        mesh = meshio.gmsh.read(str(mesh_file_name))
        meshio.vtk.write(str(vtk_file_name), mesh, binary=False, fmt_version="4.2")

        print(
            "[INFO] : VTK file for internal mesh file generated at ",
            str(mesh_file_name),
        )

    elif self.mesh_generation_method == "external":

        vtk_file_name = Path(self.output_folder) / "external.vtk"

        # Use the internal mesh to generate the vtk file
        mesh = meshio.read(str(self.mesh_file_name))
        meshio.vtk.write(str(vtk_file_name), mesh, binary=False, fmt_version="4.2")

        print(
            "[INFO] : VTK file for external mesh file generated at ",
            str(vtk_file_name),
        )

get_test_points()

This function is used to extract the test points from the given mesh

Returns:

Name Type Description
test_points ndarray

Array of test points
Source code in scirex/core/sciml/geometry/geometry_2d.py 
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
def get_test_points(self):
    """
    This function is used to extract the test points from the given mesh

    Args:
        None

    Returns:
        test_points (np.ndarray): Array of test points
    """

    if self.mesh_generation_method == "internal":
        # vtk_file_name  = Path(self.output_folder) / "internal.vtk"
        # code over written to plot from np.linspace instead of vtk file
        # generate linspace of points in x and y direction based on x and y limits
        x = np.linspace(self.x_limits[0], self.x_limits[1], self.n_test_points_x)
        y = np.linspace(self.y_limits[0], self.y_limits[1], self.n_test_points_y)
        # generate meshgrid
        x_grid, y_grid = np.meshgrid(x, y)
        # stack the points
        self.test_points = np.vstack([x_grid.flatten(), y_grid.flatten()]).T

        return self.test_points

    elif self.mesh_generation_method == "external":
        vtk_file_name = Path(self.output_folder) / "external.vtk"

    mesh = meshio.read(str(vtk_file_name))
    points = mesh.points
    return points[:, 0:2]  # return only first two columns

plot_adaptive_mesh(cells_list, area_averaged_cell_loss_list, epoch, filename='cell_residual')

Plots the residuals in each cell of the mesh.

Parameters:

Name Type Description Default
cells_list

List of cell vertices

 required
area_averaged_cell_loss_list

List of area averaged cell loss

 required
epoch

The epoch number

 required
filename

The output filename

 'cell_residual'

Returns:

Type Description

None
Source code in scirex/core/sciml/geometry/geometry_2d.py 
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
def plot_adaptive_mesh(
    self, cells_list, area_averaged_cell_loss_list, epoch, filename="cell_residual"
):
    """
    Plots the residuals in each cell of the mesh.

    Args:
        cells_list: List of cell vertices
        area_averaged_cell_loss_list: List of area averaged cell loss
        epoch: The epoch number
        filename: The output filename

    Returns:
        None
    """

    plt.figure(figsize=(6.4, 4.8), dpi=300)

    # normalise colors
    norm = mcolors.Normalize(
        vmin=np.min(area_averaged_cell_loss_list),
        vmax=np.max(area_averaged_cell_loss_list),
    )

    # Create a colormap
    colormap = plt.cm.jet

    for index, cell in enumerate(cells_list):
        x = cell[:, 0]
        y = cell[:, 1]

        x = np.append(x, x[0])
        y = np.append(y, y[0])

        curr_cell_loss = float(area_averaged_cell_loss_list[index])

        color = colormap(norm(curr_cell_loss))

        plt.fill(x, y, color=color, alpha=0.9)

        plt.plot(x, y, "k")

        # # compute x_min, x_max, y_min, y_max
        # x_min = np.min(x)
        # x_max = np.max(x)
        # y_min = np.min(y)
        # y_max = np.max(y)

        # # compute centroid of the cells
        # centroid = np.array([np.mean(x), np.mean(y)])

        # plot the loss text within the cell
        # plt.text(centroid[0], centroid[1], f"{curr_cell_loss:.3e}", fontsize=16, horizontalalignment='center', verticalalignment='center')

    sm = plt.cm.ScalarMappable(cmap=colormap, norm=norm)
    sm.set_array([])
    plt.colorbar(sm)

    # output filename
    output_filename = Path(f"{self.output_folder}/{filename}_{epoch}.png")
    plt.title(f"Cell Residual")
    plt.savefig(str(output_filename), dpi=300)

read_mesh(mesh_file, boundary_point_refinement_level, bd_sampling_method, refinement_level)

Reads mesh from a Gmsh .msh file and extracts cell information.

Parameters:

Name Type Description Default
mesh_file str

Path to the mesh file

 required
boundary_point_refinement_level int

Level of boundary point refinement

 required
bd_sampling_method str

Method for boundary point sampling ('uniform'/'lhs')

 required
refinement_level int

Level of mesh refinement

 required

Returns:

Name Type Description
cell_points

Array of cell vertices
bd_dict

Dictionary of boundary points

Raises:

Type Description
ValueError

If mesh file format is invalid
Source code in scirex/core/sciml/geometry/geometry_2d.py 
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
def read_mesh(
    self,
    mesh_file: str,
    boundary_point_refinement_level: int,
    bd_sampling_method: str,
    refinement_level: int,
):
    """
    Reads mesh from a Gmsh .msh file and extracts cell information.

    Args:
        mesh_file: Path to the mesh file
        boundary_point_refinement_level: Level of boundary point refinement
        bd_sampling_method: Method for boundary point sampling ('uniform'/'lhs')
        refinement_level: Level of mesh refinement

    Returns:
        cell_points: Array of cell vertices
        bd_dict: Dictionary of boundary points

    Raises:
        ValueError: If mesh file format is invalid
    """

    self.mesh_file_name = mesh_file

    # bd_sampling_method = "uniform"  # "uniform" or "lhs"

    file_extension = Path(mesh_file).suffix

    if file_extension != ".mesh":
        raise ValueError("Mesh file should be in .mesh format.")

    # Read mesh using meshio
    self.mesh = meshio.read(mesh_file)

    if self.mesh_type == "quadrilateral":
        # Extract cell information
        cells = self.mesh.cells_dict["quad"]

    num_cells = cells.shape[0]
    print(f"[INFO] : Number of cells = {num_cells}")
    cell_points = self.mesh.points[cells][
        :, :, 0:2
    ]  # remove the z coordinate, which is 0 for all points

    # loop over all cells and rearrange the points in anticlockwise direction
    for i in range(num_cells):
        cell = cell_points[i]
        # get the centroid of the cell
        centroid = np.mean(cell, axis=0)
        # get the angle of each point with respect to the centroid
        angles = np.arctan2(cell[:, 1] - centroid[1], cell[:, 0] - centroid[0])
        # sort the points based on the angles
        cell_points[i] = cell[np.argsort(angles)]

    # Extract number of points within each cell
    print(f"[INFO] : Number of points per cell = {cell_points.shape}")

    # Collect the Boundary point id's within the domain
    boundary_edges = self.mesh.cells_dict["line"]

    # Using the point id, collect the coordinates of the boundary points
    boundary_coordinates = self.mesh.points[boundary_edges]

    # Number of Existing Boundary points
    print(
        f"[INFO] : Number of Bound points before refinement = {np.unique(boundary_coordinates.reshape(-1,3)).shape[0] * 0.5 + 1}"
    )

    # now Get the physical tag of the boundary edges
    boundary_tags = self.mesh.cell_data["medit:ref"][0]

    # Generate a Dictionary of boundary tags and boundary coordinates
    # Keys will be the boundary tags and values will be the list of coordinates
    boundary_dict = {}

    # refine the boundary points based on the number of boundary points needed
    for i in range(boundary_coordinates.shape[0]):
        p1 = boundary_coordinates[i, 0, :]
        p2 = boundary_coordinates[i, 1, :]

        if bd_sampling_method == "uniform":
            # take the current point and next point and then perform a uniform sampling
            new_points = np.linspace(
                p1, p2, pow(2, boundary_point_refinement_level) + 1
            )
        elif bd_sampling_method == "lhs":
            # take the current point and next point and then perform a uniform sampling
            new_points = lhs(2, pow(2, boundary_point_refinement_level) + 1)
            new_points[:, 0] = new_points[:, 0] * (p2[0] - p1[0]) + p1[0]
            new_points[:, 1] = new_points[:, 1] * (p2[1] - p1[1]) + p1[1]
        else:
            print(
                f"Invalid sampling method {bd_sampling_method} in {self.__class__.__name__} from {__name__}."
            )
            raise ValueError("Sampling method should be either uniform or lhs.")

        # get the boundary tag
        tag = boundary_tags[i]

        if tag not in boundary_dict:
            boundary_dict[tag] = new_points
        else:
            current_val = new_points
            prev_val = boundary_dict[tag]
            final = np.vstack([prev_val, current_val])
            boundary_dict[tag] = final

    # get unique
    for tag in boundary_dict.keys():
        val = boundary_dict[tag]
        val = np.unique(val, axis=0)
        boundary_dict[tag] = val

    self.bd_dict = boundary_dict
    # print the new boundary points  on each boundary tag (key) in a tabular format

    total_bound_points = 0
    print(f"| {'Boundary ID':<12} | {'Number of Points':<16} |")
    print(f"| {'-'*12:<12}---{'-'*16:<16} |")
    for k, v in self.bd_dict.items():
        print(f"| {k:<12} | {v.shape[0]:<16} |")
        total_bound_points += v.shape[0]

    print(f"[INFO] : No of bound pts after refinement:  {total_bound_points}")

    # Assign to class values
    self.cell_points = cell_points

    # generate testvtk
    self.generate_vtk_for_test()

    return cell_points, self.bd_dict

write_vtk(solution, output_path, filename, data_names)

Writes the data to a VTK file.

Parameters:

Name Type Description Default
solution ndarray

The solution data to be written

 required
output_path str

The output path for the VTK file

 required
filename str

The name of the output file

 required
data_names list

List of data names

 required

Returns:

Type Description

None
Source code in scirex/core/sciml/geometry/geometry_2d.py 
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
def write_vtk(
    self, solution: np.ndarray, output_path: str, filename: str, data_names: list
):
    """
    Writes the data to a VTK file.

    Args:
        solution: The solution data to be written
        output_path: The output path for the VTK file
        filename: The name of the output file
        data_names: List of data names

    Returns:
        None
    """
    # read the existing vtk into file
    if self.mesh_generation_method == "internal":
        vtk_file_name = Path(self.output_folder) / "internal.vtk"
    elif self.mesh_generation_method == "external":
        vtk_file_name = Path(self.output_folder) / "external.vtk"

    data = []
    with open(vtk_file_name, "r", encoding="utf-8") as File:
        for line in File:
            data.append(line)

    # get the output file name
    output_file_name = Path(output_path) / filename

    if solution.shape[1] != len(data_names):
        print("[Error] : File : geometry_2d.py, Function: write_vtk")
        print(
            "Num Columns in solution = ",
            solution.shape[1],
            " Num of data names = ",
            len(data_names),
        )
        raise ValueError("Number of data names and solution columns are not equal")

    # write the data to the output file
    with open(str(output_file_name), "w", encoding="utf-8") as FN:
        for line in data:
            FN.write(line)
            if "POINT_DATA" in line.strip():
                break

        for i in range(solution.shape[1]):
            FN.write("SCALARS " + data_names[i] + " float\n")
            FN.write("LOOKUP_TABLE default\n")
            np.savetxt(FN, solution[:, i])
            FN.write("\n")

print_table(title, columns, col_1_values, col_2_values)

This function prints a table with two columns to the console.

Parameters:

Name Type Description Default
title str

str: Title of the table

 required
columns list

list: List of column names

 required
col_1_values list

list: List of values for column 1

 required
col_2_values list

list: List of values for column

 required

Returns:

Type Description

None
Source code in scirex/core/sciml/utils/print_utils.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
def print_table(title: str, columns: list, col_1_values: list, col_2_values: list):
    """
    This function prints a table with two columns to the console.

    Args:
        title: str: Title of the table
        columns: list: List of column names
        col_1_values: list: List of values for column 1
        col_2_values: list: List of values for column

    Returns:
        None
    """

    # Create a console object
    console = Console()

    # Create a table with a title
    table = Table(show_header=True, header_style="bold magenta", title=title)

    # Add columns to the table
    for column in columns:
        table.add_column(column)

    # Add rows to the table
    for val_1, val_2 in zip(col_1_values, col_2_values):
        # check if val_2 is a float
        if isinstance(val_2, float):
            # add the row to the table
            table.add_row(val_1, f"{val_2:.4f}")
        else:
            table.add_row(val_1, str(val_2))

    # Print the table to the console
    console.print(table)