Hello all,
I am unable to edit my original post. So I am creating a new post with the
equations. I think that I solved issue with the equations not being
displayed. The original posting will be deleted.
I have recently started to learn how to use Deal.II. I have been spending
some going through some of the video tutorials for Deal.II, going through
the coding tutorials, and watching some of the videos from the University
of Michigan related to Finite Element Analysis. I have decided that I need
to get my hands a bit dirty so that everything that I am reading and
watching will start to sink.
So please, bear with me. I am still quite new to this field and I am still
learning some of the basics.
With that being said, I am looking at creating my own electrostatic solver.
For those unfamiliar, I am solving the PDEs:
For the electrostatic simulation, I would be solving for the Voltage scalar
field and the Electric field (which is a vector).
So far, I have found that there are two paths that I can take to get my
bilinear form. I have already attempted one path and I ran into a compile
error with Deal.II. So, I am thinking that this bilinear form may not be
appropriate. But I wanted to turn to the community here to see if my
bilinear form has an issue or if my code has an issue.
For the first path, since V and E are my unknowns, I wanted both V and E to
be on the left-hand side of their PDEs. With a bit of substitution, the
equations that I would be solving turned into this:
Normally, there would be a - signed on the right-hand side of both
equations. This will get dropped out because I am going to simplify the
equation a bit further later down the road. In order to construct the
bilinear form, I used the technique introduced in lecture 19. I am not
going to show the entire process but rather the end result. If the entire
process is needed, please post and I will respond with the complete
step-by-step process that I performed.
then
This is how I created my bilinear form. Note that the p/e term dropped out
because in this simplified simulation, there is no charge density. I create
my simulation that is based on the tutorial step-20. There are alot of
extra features in the tutorial that I did not add in to my simulation but
rather focused on the FESystem portion. My solver is a CG Solver. Nothing
fancy here. My first version of the solver can be found below
void SetupSolver(FESystem<2> &finiteElement, DoFHandler<2> &dofHandler,
SparseMatrix<double> &masterMatrix, Vector<double> &masterRHS, Vector
<double> &masterSolution)
{
QGauss<2> quadrature(3);// Sets the form of the function
FEValues<2> finiteElementValues(finiteElement, quadrature,
update_values | update_gradients | update_JxW_values |
update_quadrature_points);
const unsigned int dofsPerCell = finiteElement.dofs_per_cell;
const unsigned int numberQuadraturePoints = quadrature.size();
FullMatrix<double> cell_matrix(dofsPerCell, dofsPerCell);
Vector<double> rhsCell(dofsPerCell);
vector<types::global_dof_index> localDofIndices(dofsPerCell);
DoFHandler<2>::active_cell_iterator cell = dofHandler.begin_active();
DoFHandler<2>::active_cell_iterator endc = dofHandler.end();
const FEValuesExtractors::Vector electricField(2);
const FEValuesExtractors::Scalar voltageField(0);
cout << "Setting up Solver" << endl;
for(; cell != endc; cell++)
{
// Calculates the values and gradients of the shape function for a
particular cell
finiteElementValues.reinit(cell);
cell_matrix = 0;
rhsCell = 0;
cout << "working on cell: " << cell << endl;
for(unsigned int qIndex = 0; qIndex < numberQuadraturePoints; qIndex
++)
{
for(unsigned int i = 0; i < dofsPerCell; i++)
{
const double div_phi_i_e = finiteElementValues[electricField
].divergence(i, qIndex);
const Tensor<1, 2> grad_phi_i_v = finiteElementValues[
voltageField].gradient(i, qIndex);
for(unsigned int j = 0; j < dofsPerCell; j++)
{
const Tensor<1, 2> phi_j_e = finiteElementValues[
electricField].value(j, qIndex);
const Tensor<1, 2> grad_phi_j_v = finiteElementValues[
voltageField].gradient(j, qIndex);
cell_matrix(i, j) += ((div_phi_i_e * phi_j_e) +
(grad_phi_i_v
* grad_phi_j_v)) * finiteElementValues.JxW(qIndex);
}
rhsCell(i) += 0.0;
}
}
cell->get_dof_indices(localDofIndices);
for(unsigned int i = 0; i < dofsPerCell; i++)
{
for(unsigned int j = 0; j < dofsPerCell; j++)
{
masterMatrix.add(localDofIndices[i], localDofIndices[j],
cell_matrix(i, j));
}
}
for(unsigned int i = 0; i < dofsPerCell; i++)
{
masterRHS(localDofIndices[i]) += rhsCell(i);
}
}
cout << "Finished setting up cells" << endl;
// SolvingFunction test();
map<types::global_dof_index, double> boundaryValues;
cout << "Interperlating boundary conditions" << endl;
VectorTools::interpolate_boundary_values(dofHandler, 1, ConstantFunction
<2>(100.), boundaryValues);
VectorTools::interpolate_boundary_values(dofHandler, 2, ZeroFunction<2
>(), boundaryValues);
VectorTools::interpolate_boundary_values(dofHandler, 0, ZeroFunction<2
>(), boundaryValues);
MatrixTools::apply_boundary_values(boundaryValues, masterMatrix,
masterSolution, masterRHS);
cout << "Finished setting up solver" << endl;
}
On line
cell_matrix(i, j) += ((div_phi_i_e * phi_j_e) + (grad_phi_i_v * grad_phi_j_v
)) * finiteElementValues.JxW(qIndex);
I am getting the error:
error: no match for 'operator+' (operand types are 'dealii::Tensor<1, 2>'
and 'dealii::Tensor<0, 2, double>::tensor_type {aka double}')
Which, from what the error is stating, the + operator can not take these
two arguments.
This then leads me to reason that my bilinear form has an issue. Would this
be a correct assumption? Or did I code something incorrectly? (This is my
main issue right now)
Path 2:
I have also been considering this approach for creating the bilinear form
where the PDEs that I want to solve are as follows:
Again, following the technique introduced in Lecture 19, I would simplify
my PDEs into:
then
I was wondering if this would be a better approach in constructing the
bi-linear form and solving for the electric field and voltage field.
Or are both methods equally valid? I have a feeling that the first method
is not valid since it is essentially 2 of the same equations. What are your
thoughts? Is there a different method that I can approach in solving the
PDEs for the electric and voltage fields?
Side Notes:
There is a third equation which is the curl of the electric field is 0 for
electrostatics. I decided that I am not using it since these two are
enough. Is this an incorrect thinking and should I also account for this
equation in my bilinear form?
I am using Deal.II to create a mesh which is a parallel plate capacitor. I
am then solving for the electric and voltage fields that make up the
parallel plates when the voltage and electric fields are 0 on the boundary.
On the boundary of one of the plates, the Voltage is set to 100. On the
other plate, the Voltage is set to 0.
For reference, I am attaching the source code of my simulation. There are
more improvements that I will be making to the code once I get something
workable. The code that is attached is a modified version of an earlier
simulation where I was only solving for the voltage field in a parallel
plate. This was around 2 years ago that I created this. I am simply
expanding the code to also solve for the electric field. in short, it is
rough. If you do see any other mistakes or areas of improvement, feel free
to let me know.
Thank you
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/*
This file is the previous version of the capacitor simulator.
This builds and does simulate!
*/
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_values_extractors.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <fstream>
#include <cmath>
#include <set>
using namespace dealii;
using namespace std;
using namespace GridGenerator;
// Generates the grid needed
void generateGrid(unsigned int radius, double distance, double height, double length, Triangulation<2> &finalProblem)
{
const Point<2> centerBoundary (0,0);
const SphericalManifold<2> manifold_description(centerBoundary);
set<Triangulation<2>::active_cell_iterator> cellsToDelete;
Triangulation<2>::active_cell_iterator cell;
Triangulation<2>::active_cell_iterator endCell;
GridOut outputGrid;
// This is the Triangulation of the entire problem
Triangulation<2> parallelPlate;
hyper_ball(parallelPlate, centerBoundary, radius);
cout << "Creating Grid" << endl;
if(length > distance)
{
cout << "This will NOT work";
return;
}
parallelPlate.set_all_manifold_ids(0);
for (auto cell : parallelPlate.active_cell_iterators())
{
if (cell->center().square() < 1.e-10)
{
cell->set_manifold_id(numbers::invalid_manifold_id);
for(unsigned int i = 0; i < GeometryInfo<2>::faces_per_cell; ++i)
{
cell->face(i)->set_manifold_id(numbers::invalid_manifold_id);
}
}
}
cell = parallelPlate.begin_active();
endCell = parallelPlate.end();
parallelPlate.set_manifold(0, manifold_description);
parallelPlate.refine_global(4);
//TODO: add in flag for smoothing-> Triangulation::MeshSmoothing
// This will interate over all of the cells and delete the cells to make the two parallel plates
cout << "Creating parallel Plates" << endl;
// Iterates through all of the cells
for(;cell != endCell; cell++)
{
for(unsigned int v = 0; v < 4; v++)
{
if(abs(cell->vertex(v)(1)) <= (height / 2) && abs(cell->vertex(v)(0)) >= (distance / 2) && abs(cell->vertex(v)(0)) <= (distance / 2 + length))
{
// Marks the cells that fit the if statement for deletion
cellsToDelete.insert(cell);
break;
}
}
}
parallelPlate.set_manifold(0);
GridGenerator::create_triangulation_with_removed_cells(parallelPlate, cellsToDelete, finalProblem);
// Reset the cell and endCell interators
cell = finalProblem.begin_active();
endCell = finalProblem.end();
// This will interrate over all of the cells and refine the ones that are within the boundaries of the plates within the air
cout << "Refining grid around plates" << endl;
for(int step = 0; step < 2; step++)
{
// Iterates through all of the cells
for(;cell != endCell; cell++)
{
for(unsigned int v = 0; v < 2; v++)
{
if(abs(cell->vertex(v)(1)) <= 4.4 && abs(cell->vertex(v)(0)) <= 4.4)
{
// Marks the cells that fit the if statement for refinement
cell->set_refine_flag();
break;
}
}
}
cell = finalProblem.begin_active();// Resets the cell back to the beginning
// Creates a new mesh with the refined cells
finalProblem.execute_coarsening_and_refinement();
}
cell = finalProblem.begin_active();
endCell = finalProblem.end();
// Iterate through all of the cells and determine which ones are at the boundary for the plates
cout << "Determining Boundary of parallel Plates" << endl;
for(; cell != endCell; cell++)
{
if(cell->at_boundary())
{
for(unsigned int f = 0; f < GeometryInfo<2>::faces_per_cell; f++)
{
if(cell->face(f)->at_boundary())
{
for(unsigned int v = 0; v < GeometryInfo<2>::vertices_per_face; v++)
{
// Change the boundary ID for the cell here
// These two checks will ensure that only the cells related to the edge of the boundary for the parallel plates will be selected
// and modified
// The + 1 is added to ensure that the cells on the back side of the boundary are selected
if(abs(cell->vertex(v)(0)) <= (distance / 2 + length + 1) && abs(cell->vertex(v)(1)) <= (height) && cell->vertex(v)(0) < 0)
{
// Sets the Boundary ID for the left plate
cell->face(f)->set_boundary_id(1);
}
else if(abs(cell->vertex(v)(0)) <= (distance / 2 + length + 1) && abs(cell->vertex(v)(1)) <= height && cell->vertex(v)(0) > 0)
{
// Sets the Boundary ID for Right plate
cell->face(f)->set_boundary_id(2);
}
}
}
}
}
}
ofstream out ("plate.vtk");
outputGrid.write_vtk(finalProblem, out);
cout << "Finish creating grid" << endl;
}
void GenerateDOFS(DoFHandler<2> &dofHandler, FESystem<2> &finiteElement, SparsityPattern &sparsePattern, SparseMatrix<double> &masterMatrix, Vector<double> &masterRHS, Vector<double> &masterSolution)
{
dofHandler.distribute_dofs(finiteElement);
cout << "calculating degrees of freedom" << endl;
// Creates a dynamic Sparsity Pattern so that the sparsity pattern function does not over estimate the number of elements
DynamicSparsityPattern dynamicSparsityPattern(dofHandler.n_dofs(), dofHandler.n_dofs());
DoFTools::make_sparsity_pattern(dofHandler, dynamicSparsityPattern);
sparsePattern.copy_from(dynamicSparsityPattern);// Copies the dynamic sparsity pattern into the sparse pattern object
cout << "Number of degress of freedom: " << dofHandler.n_dofs() << endl;
masterMatrix.reinit(sparsePattern);// Initilizes the masterMatrix with the sparse pattern
masterRHS.reinit(dofHandler.n_dofs());
masterSolution.reinit(dofHandler.n_dofs());
}
void SetupSolver(FESystem<2> &finiteElement, DoFHandler<2> &dofHandler, SparseMatrix<double> &masterMatrix, Vector<double> &masterRHS, Vector<double> &masterSolution)
{
QGauss<2> quadrature(3);// Sets the form of the function
FEValues<2> finiteElementValues(finiteElement, quadrature, update_values | update_gradients | update_JxW_values | update_quadrature_points);
const unsigned int dofsPerCell = finiteElement.dofs_per_cell;
const unsigned int numberQuadraturePoints = quadrature.size();
FullMatrix<double> cell_matrix(dofsPerCell, dofsPerCell);
Vector<double> rhsCell(dofsPerCell);
vector<types::global_dof_index> localDofIndices(dofsPerCell);
DoFHandler<2>::active_cell_iterator cell = dofHandler.begin_active();
DoFHandler<2>::active_cell_iterator endc = dofHandler.end();
const FEValuesExtractors::Vector electricField(2);
const FEValuesExtractors::Scalar voltageField(0);
cout << "Setting up Solver" << endl;
for(; cell != endc; cell++)
{
// Calculates the values and gradients of the shape function for a particular cell
finiteElementValues.reinit(cell);
cell_matrix = 0;
rhsCell = 0;
cout << "working on cell: " << cell << endl;
for(unsigned int qIndex = 0; qIndex < numberQuadraturePoints; qIndex++)
{
for(unsigned int i = 0; i < dofsPerCell; i++)
{
const double div_phi_i_e = finiteElementValues[electricField].divergence(i, qIndex);
const Tensor<1, 2> grad_phi_i_v = finiteElementValues[voltageField].gradient(i, qIndex);
for(unsigned int j = 0; j < dofsPerCell; j++)
{
const Tensor<1, 2> phi_j_e = finiteElementValues[electricField].value(j, qIndex);
const Tensor<1, 2> grad_phi_j_v = finiteElementValues[voltageField].gradient(j, qIndex);
cell_matrix(i, j) += ((div_phi_i_e * phi_j_e) + (grad_phi_i_v * grad_phi_j_v)) * finiteElementValues.JxW(qIndex);
}
rhsCell(i) += 0.0;
}
}
cell->get_dof_indices(localDofIndices);
for(unsigned int i = 0; i < dofsPerCell; i++)
{
for(unsigned int j = 0; j < dofsPerCell; j++)
{
masterMatrix.add(localDofIndices[i], localDofIndices[j], cell_matrix(i, j));
}
}
for(unsigned int i = 0; i < dofsPerCell; i++)
{
masterRHS(localDofIndices[i]) += rhsCell(i);
}
}
cout << "Finished setting up cells" << endl;
// SolvingFunction test();
map<types::global_dof_index, double> boundaryValues;
cout << "Interperlating boundary conditions" << endl;
VectorTools::interpolate_boundary_values(dofHandler, 1, ConstantFunction<2>(100.), boundaryValues);
VectorTools::interpolate_boundary_values(dofHandler, 2, ZeroFunction<2>(), boundaryValues);
VectorTools::interpolate_boundary_values(dofHandler, 0, ZeroFunction<2>(), boundaryValues);
MatrixTools::apply_boundary_values(boundaryValues, masterMatrix, masterSolution, masterRHS);
cout << "Finished setting up solver" << endl;
}
void solve(SparseMatrix<double> &masterMatrix, Vector<double> &masterSolution, Vector<double> &masterRHS)
{
SolverControl solverControl(100000, 1e-12);
SolverCG<> solver(solverControl);
cout << "Beginning solver" << endl;
solver.solve(masterMatrix, masterSolution, masterRHS, PreconditionIdentity());
cout << "Finished solving" << endl;
}
void outputResults(DoFHandler<2> &masterDOFHandler, Vector<double> &masterSolution)
{
std::vector<std::string> solutionNames(2, "Electric Field");
solutionNames.push_back("Voltage Value");
std::vector<DataComponentInterpretation::DataComponentInterpretation> interpret(2, DataComponentInterpretation::component_is_part_of_vector);
interpret.push_back(DataComponentInterpretation::component_is_scalar);
DataOut<2> outputData;
outputData.add_data_vector(masterDOFHandler, masterSolution, solutionNames, interpret);
outputData.build_patches();
ofstream output("solution.vtk");
outputData.write_vtk(output);
}
int main()
{
Triangulation<2> testtri;
FE_Q<2> finiteElementQuadrature(2);
FESystem<2> electrostaticSystem;
SparsityPattern sparsePattern;
SparseMatrix<double> systemMatrix;
Vector<double> solution;
Vector<double> RHS;
generateGrid(15, 2, 5, 0.5, testtri);
DoFHandler<2> internalDOFHandler(testtri);
GenerateDOFS(internalDOFHandler, electrostaticSystem, sparsePattern, systemMatrix, RHS, solution);
SetupSolver(electrostaticSystem, internalDOFHandler, systemMatrix, RHS, solution);
solve(systemMatrix, solution, RHS);
outputResults(internalDOFHandler, solution);
return 0;
}
