|
last change: 04.03.2008 New solver - Solving methods The main method for solving, which is called by the NeuroFEM package, is the Solve method. The task of this method is to acquire the right hand side for the problem, to initialize the solution, to solve the problem, to store the result back to NeuroFEM and to asure the value at the reference electrode position. For each of those parts a seperate method is used, which will all be described below (except the SolveProblem method which can be found on the next page).Code from file "neurofem/fmg/fmgsolver.cpp" void FMGSolver::Solve() {InitValues(); SolveProblem(); StoreSolution(); ShiftSolution(); } The InitValues method initializes the solution and loads the right hand side of the problem from NeuroFEM. For this task some of the already described methods can be used (especially the mapping methods). Code from file "neurofem/fmg/fmgsolver.cpp" void FMGSolver::InitValues() {// Initialize the values *u = 0.0; *f = 0.0; // Copy right-hand-side from neurofem int i,j,k; ParExPDE::CuboidStorageBase for (int l = 1; l <= numUnknowns; l++) { MapNodeToIndices(l, &i,&j,&k); f_array->valueAt(i,j,k) = -rhs[l-1]; } } The StoreSolution method stores the informations from ParExPDE back to NeuroFEM. For this the same mapping methods are used as for the InitValues method. Code from file "neurofem/fmg/fmgsolver.cpp" void FMGSolver::StoreSolution() {// Copy the solution to neurofem int i,j,k; ParExPDE::CuboidStorageBase for (int l = 1; l < numUnknowns; l++) { MapNodeToIndices(l, &i,&j,&k); sol[l-1] = u_array->valueAt(i,j,k); } } The ShiftSolution method asures the correctness of the solution value at the reference electrode. For this task the obtained solution at the position is compared with the value of the reference electrode and the whole solution is shifted to match both values. This shifting is allowed due to the pure Neumann boundary conditions. Code from file "neurofem/fmg/fmgsolver.cpp" void FMGSolver::ShiftSolution() {int* dir_index = VolumeGrid.getGridNodeDirichletIndex(); double* dir_value = VolumeGrid.getGridNodeDirichletPotential(); int index = -1; double value = 0.0; for (int i = 0; i < numUnknowns; i++) if (dir_index[i] > 0) { index = i; value = dir_value[i]; } // Do the shifting if (index != -1) { double shift = sol[index] - value; for (int i = 0; i < numUnknowns; i++) sol[i] -= shift; } } | |