hello, I am writing topology optimization code and want to collect data using DataFrame but the result as shown in DataFrame is only the data in the last iteration , how can DataFrame collect the result in every iteration.
import dolfinx
import numpy as np
import pandas as pd
import dolfinx.nls.petsc
import ufl
import dolfinx.fem.petsc
import pyvista
import matplotlib.pyplot as plt
from dolfinx.mesh import CellType, create_rectangle
from mpi4py import MPI
from dolfinx import fem, nls, io, default_scalar_type
from ufl import Measure, inner, grad, div, TestFunction, TrialFunction, adjoint, derivative, SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx
from dolfinx.fem import (Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from petsc4py import PETSc
from dolfinx.plot import vtk_mesh
import pyvista as pv
import time
#Library for ML
import seaborn as sns
import tensorflow as tf
import tensorflow.compat.v1 as tf
from sklearn.model_selection import train_test_split
from tensorflow.keras.layers import Dense
from tensorflow.keras.models import Sequential
from sklearn.preprocessing import LabelEncoder
import ast
start_time = time.time()
#Parameter
D0 = 1;
k0 = -1.0;
eta = 2.0;
gamma = 2.0;
alpha = 1e-16
Rfl = 0.0018
# mesh and function space
mesh = create_rectangle(MPI.COMM_WORLD, [[0, 0], [1, 1]], [100, 100], CellType.quadrilateral)
V = fem.FunctionSpace(mesh, ("CG", 1))
c = fem.Function(V)
phi = TestFunction(V)
lambdh = fem.Function(V)
lamda_real,sen_real,eps_real,c_real = fem.Function(V),fem.Function(V),fem.Function(V),fem.Function(V)
# The eps is the design variable
eps = fem.Function(V)
eps.x.array[:] = 0.5
data_sen = []
data_lam = []
data_eps = []
data_c = []
data_ite = []
loop = 5
for i in range(loop):
eps_real = eps.x.array.real
data_eps.append(eps_real)
df_eps = pd.DataFrame(data_eps)
df_eps_final = pd.DataFrame(df_eps)
print("iteration ",i)
print("eps",eps_real)
#boundary_markers = FacetFunction('size_t',mesh)
ds = Measure('ds')
dx = Measure("dx") ##integration over different parts of the domain
x = SpatialCoordinate(mesh)
gu = Constant(mesh, default_scalar_type(0)) ##Neauman = 0
# Residual of the variational form of Poisson problem
R = inner(D0*eps**eta*grad(c), grad(phi))*dx - inner(k0*(1-eps)**gamma*c,phi)*dx - inner(gu,phi) *ds # " - inner(gu,phi) *ds " is boundary neauman
##Dirichlet
def on_boundary(x):
return np.isclose(x[0], 0)
dofs_L = fem.locate_dofs_geometrical(V , on_boundary)
bc_L = fem.dirichletbc(default_scalar_type(1), dofs_L, V) ##Dirichlet = 1
bcs = [bc_L]
#Forming and solving the linear system
problem_u = fem.petsc.NonlinearProblem(R, c, bcs)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD,problem_u)
solver.rtol = 1e-16
solver.solve(c)
#how big is the solution
J = -ufl.inner((k0*(1-eps)**gamma),c)*dx ##Objective function
J_form = fem.form(J)
J_old = fem.assemble_scalar(J_form)
#solve the adjoint vector
lhs = ufl.adjoint(ufl.derivative(R,c))
rhs = -ufl.derivative(J,c)
def on_boundary_1(x):
return np.isclose(x[0], 0)
dofs_L_adj = fem.locate_dofs_geometrical(V , on_boundary_1)
bc_L_adj = fem.dirichletbc(default_scalar_type(0), dofs_L_adj, V) ##Dirichlet = 1
bcs_adjoint = [bc_L_adj]
#Adjective problem
problem_adj = dolfinx.fem.petsc.LinearProblem(lhs, rhs, bcs_adjoint, petsc_options={"ksp_type": "preonly", "pc_type": "lu","pc_factor_mat_solver_type": "mumps"})
#solve adj problem in lamda (lamda is adjoint variable)
lamda = problem_adj.solve()
dJde = dolfinx.fem.Expression(k0*gamma*((1-eps)**(gamma-1))*c+dot(D0*eta*(eps**(eta-1))*grad(c),grad(lamda))+inner(k0*gamma*((1-eps)**(gamma-1))*c,lamda), V.element.interpolation_points())
sens = fem.Function(V)
sens.interpolate(dJde)
# steepest descent with constant beta as stepsize
beta = 0.02
delta_eps = beta*sens
eps_new = dolfinx.fem.Expression(eps + delta_eps, V.element.interpolation_points())
# determine new value for the design variables
eps.interpolate(eps_new)
eps1 = fem.Function(V)
phi_3 = TestFunction(V)
filter = inner((Rfl**2)*grad(eps1), grad(phi_3))*dx + inner(eps1,phi_3)*dx - inner(eps,phi_3)*dx
filter_eps = fem.petsc.NonlinearProblem(filter, eps1)
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD,filter_eps)
solver.rtol = 1e-16
solver.solve(eps1)
eps.interpolate(eps1)
eps.x.array[eps.x.array > 1] = 1
eps.x.array[eps.x.array < 0] = 0
import math
lamda_real,sen_real,eps_real,c_real = lamda.x.array.real, sens.x.array.real, eps.x.array.real, c.x.array.real
###collect data###
data_c.append(c_real)
df_c = pd.DataFrame(data_c)
df_c_final = pd.DataFrame(df_c)
print("c",c_real)
print(" ")
#data = eps
#df = pd.DataFrame(data)
end_time = time.time()
execution_time = end_time - start_time
print("Execution time:", execution_time, "seconds")
the result from printed c and eps
iteration 0
eps [0.5 0.5 0.5 ... 0.5 0.5 0.5]
c [1. 1. 0.99243388 ... 0.64808462 0.64805222 0.64805222]
iteration 1
eps [0.49144952 0.49144952 0.49145228 ... 0.49159995 0.4916007 0.4916007 ]
c [1. 1. 0.99202673 ... 0.63116817 0.63113441 0.63113441]
iteration 2
eps [0.4836042 0.4836042 0.48359976 ... 0.48349888 0.48350033 0.48350033]
c [1. 1. 0.99163397 ... 0.61468617 0.61465109 0.61465109]
iteration 3
eps [0.47648518 0.47648518 0.47646429 ... 0.47569298 0.47569508 0.47569508]
c [1. 1. 0.99126011 ... 0.59872791 0.59869154 0.59869154]
iteration 4
eps [0.47009843 0.47009843 0.47005283 ... 0.46817523 0.46817795 0.46817795]
c [1. 1. 0.990909 ... 0.58336524 0.5833276 0.5833276 ]
show dataframe collecting ‘c’
show dataframe collecting ‘eps’
